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CHAPTER 1 


INTRODUCTION 


The analysis of turbulent compressible jets is a. critical technology for improving 
the performance of advanced aerospace propulsion systems. In bot h the military and 
commercial sectors a need exists for a better understanding and prediction of jet 
engine exhaust plumes and fuel injector flows. This dissertation addresses this need 
through the development and application of a computational method to analyze these 
flows, providing a better understanding of the underlying physics. 

1.1 Motivation 

The survivability of military aircraft depends on their ability to evade the enemy s 
defenses. Anti-aircraft systems employ infrared detectors and heat, seeking missiles, 
which rely on the high temperature exhaust from the aircraft’s engines. The military 
has a great interest in low-observable technology to reduce the size and intensity of 
this exhaust plume in order to increase the aircraft's survival rate. 

In the commercial aircraft industry, reducing the noise generated by the aircraft 
has become an important focus, as increasingly stringent restrictions on noise levels 
are imposed near airports. The primary contribution to community noise is jet noise 
from the engine at takeoff. The mechanisms by which the jet generates sound are 
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nol well understood. Expanding the knowledge of jet aeroaeoust ics lias become the 
focus of several recent national programs including NASA's High Speed Research and 
Advanced Subsonic Technology programs. 

Both the military arid commercial sectors are interested in high speed ramjet and 
sc ramjet- propulsion systems. These systems are being considered for fast response 
weapon systems, high speed aircraft and efficient airbreathing launch vehicles. The 
complete and rapid mixing of the fuel jet with the high speed airstrcam is key to 
their operation. A better predictor of the performance of fuel injectors and their 
interaction with the freestream fiowfield is necessary. 

1.2 Physics of Jets 

The physics of jet flows is dominated by turbulent motion. Turbulent flows bv 
definition are unsteady and randomly varying. Turbulence is three-dimensional and 
rotational with many vortical structures. The scales of the structures vary from the 
Kolmogorov microscales [1] to scales nearly on the order of the jet diameter. The 
large scales contain most of the turbulent energy and transport the majority of the 
momentum and energy. The energy is cascaded from the largest scales to the smallest. 
The smallest scales then dissipate the energy and are isotropic. The Kolmogorov 
scales, k . are extremely small. Their ratio to the scales of the largest eddies is 

J = fler 3/ " (U) 

where f. is the size of the largest scales and Ret is the turbulent Reynolds number. 
However, Wilcox shows for a typical turbulent flow the Kolmogorov scales are ap- 
proximately seventy two times the mean free path of the molecules [2]. Thus, the 
continuum assumption can be used when modeling the flow at this level. 
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The turbulent motion in the jet flow-field is confined to the mixing layer, which 
is the interface between the flow from the nozzle and the ambient air (figure 1.1). 
The mixing layer begins at the jet lip and spreads radially with increasing distance 
downstream. The nozzle flow that has not been affected by the mixing layer maint ains 
an inviscid character and is termed the potential core. 

1.3 Computational Methods 

An accurate and inexpensive method to predict the complex physics of jet flows 
would be a great benefit to low-observable, aeroacoustic. and fuel injector technology. 
Experimental studies of jets are an expensive and difficult undertaking and provide 
a limited amount of data. Analytical methods in this area are very limited in their 
applicability due to the steady-state and low Reynolds number assumptions that 
usually accompany them. 

Computational fluid dynamics (CFD) offers an excellent alternative for jet analy- 
ses. In CFD. solutions to the governing equations of fluid motion are obt ained using 
numerical methods. In principle, it does not require major simplifications to the 
equations that limit the applicability and accuracy of the simpler analytical methods. 
CFD provides a complete description of the flowfield at specified discretized points. 
In this way it is superior to experimental methods that provide a. limited amount of 
data.. In general, a CFD analysis is also less expensive than an experimental program. 

Two sources of error limit the accuracy of CFD. The first error is the error in- 
troduced into the solution by the discretization of the equations. This discretization, 
or truncation, error is a function of both the numerical scheme used to solve the 
equations and the computational grid that specifies the discrete locations at which 
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tilt' equations are solved. The second source of error is modeling error. This error 
is the result of assumptions made to model the particular geometry or phvsics of 
the problem. Modeling errors are introduced in grid generation, boundary condition 
specification and within problem formulation and solution process. In fact some of 
the largest and most widespread sources of modeling error occur in formulating and 
solving a turbulent flowfield. Three common methods are used to simulate turbulent 
flows. They are outlined below. 

1.3.1 Reynolds Averaged Navier-Stokes Simulations 

The most common approach to simulate a turbulent flow is to solve the Reynolds 
Averaged Navier-Stokes (RANS) equations. The RANS equations are obtained by 
time averaging the Navier-Stokes equations. The contribution of the unsteady terms 
in the equations is averaged out and the effect on the flow is replaced by the Reynolds 
stress tensor. In practice, this tensor is modeled using the Boussinesq approximation: 
it is replaced by the product of an eddy viscosity and the strain rate tensor. To 
further simplify, isotropic turbulence is typically assumed. The process of calculating 
the eddy viscosity is commonly referred to as turbulence modeling. 

Turbulence modeling has been an active area of research for many' years and a great 
number approaches exists. But due to the approximations inherent in the method, 
these approaches have failed to produce an adequate simulation of a turbulent jet. 

The most widely used turbulence model for RANS simulations of jet flows is the k 
< model [3]. In this approach, two additional partial differential equations, transport 
equations for the turbulent kinetic energy, k, and the turbulent dissipation rate, e, 
are solved and the eddy viscosity' is computed from these quantities. 
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Corrections to the k-c model to improve jet flow predictions have been developed. 
Sarkar [4] and Zemaii [5] both introduced corrections to account for compressibility 
effects which retard the shear layer mixing in high speed jets. Pope [(>] developed a 
vortex stretching correction that, improved ihe predict ion of round jets. 

In a. cooperative effort among the aerospace community. Barber et. al. [7] com- 
puted several jet flows using RANS techniques. They found that the solutions were 
highly dependent on the formulation of the turbulence model and the corrections 
used. They concluded that no one model provided adequate predictions over a range 
of jet conditions. 

RANS simulations are the least CPU intensive method for computing a turbulent 
flow. Typically, the discretized equations are marched in time to convergence at 
a steady state. Time accuracy in intermediate steps is not necessary and is often 
sacrificed for computational speed. In general, the results obtained represent only 
the time average of the flowfield. Some unsteady information may be available from 
the turbulence model itself (i.e. turbulent kinetic energy is computed in the k-c 
model). 

1.3.2 Direct Numerical Simulations 

In theory the simplest and most straightforward way to compute a turbulent 
flowfield is by performing a Direct Numerical Simulation (DNS). DNS methods solve 
the Navier-Stokes equations in a time accurate manner without approximation. In 
order to obtain an accurate representation of the turbulence, the turbulent motion 
down to the Kolmogorov scales must be accurately resolved in both time and space. 
In other words, the grid spacing must be no larger than the Kolmogorov scale, k. The 
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Kolmogorov scales vary with the turbulent Reynolds number as shown in equation 
(1.1). In order to resolve just one large scale eddy we would need k grid points in 
each direction, a total of /?<y points. And since the time step is linearly related to 
the grid size through the Courant Friedrichs Lewey (OFL) number [8]. the cost of the 
simulation is on the order of Ref. The amount of computer memory available limits 
the size of the computational grid that can be used and the computer’s processing 
speed gives a practical limit to the number of time steps possible in a simulation. 
Therefore, based on available computational resources. DNS is limited to very low 
Reynolds number flows. 

One way to help alleviate this Reynolds number limitation on DNS is the use 
of high order numerical methods [9 15]. These methods have the ability to accu- 
rately capture small scale structures with fewer grid points than possible with tradi- 
tional second-order accurate codes. This increased accuracy does come at the price 
of increased computational expense and the trade-off between the two has not been 
sufficiently investigated. 

Freund [16-18] has computed Mach 0.8. 0.9, and 1.92 jets using a DNS technique. 
These jets had Reynolds numbers of 800, 3600, and 2000 respectively. This work has 
produced excellent agreement with experimental data and provided great insight into 
the flow behavior. But the Reynolds numbers of the jets are orders of magnitude 
smaller than any nozzle of practical interest. 

1.3.3 Large-Eddy Simulations 

A compromise between the approximations necessary in RANS and the computa- 
tional limitations of DNS is Large-Eddy Simulation (LES) [19-23]. In LES the large 
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scales of turbulent motion are simulated directly in the Navier-Stokes equations, but 
the small scales are modeled in a manner similar to RANS turbulence modeling. 
Because the larger scales carry most of the momentum and energy, computing them 
directly should increase the simulation's accuracy. And. since the small scales are dis- 
sipative and isotropic, modeling t hem using a simple eddy viscosity approach appears 
viable. The scales are separated by spatially filtering the Navier-Stokes equations. 
This filtering process replaces the equations with a set of resolved (large scale) equa- 
tions of motion that contain additional unresolved (small scale) terms that must be 
modeled. The size of the scales that are resolved and that are modeled is determined 
by the width of the filter. A. which is on the order of the grid cell size. 

The key to an accurate LES computation is the model used to approximate the 
unresolved or sub-grid scale terms. There are many forms of these sub-grid scale 
models. The simplest and most popular model, the Smagorinsky model [24] is similar 
to Prandtl’s mixing length theory [25]. Germanoet. al. [26] developed a very success- 
ful model which dynamically adjusts this constant based on the local flow conditions. 
Others [27-30] have modified and improved on these basic ideas. 

As with DNS calculations, high order schemes are important for proper and effi- 
cient resolution of the large scale structures. In addition, most sub-grid models scale 
the eddy viscosity by the square of the filter width, which is directly related to the 
grid size. In second-order accurate numerical schemes the truncation error is also a 
function of the square of the cell size. LES performed with a second-order accurate 
code would have a truncation error and a sub-grid model of the same magnitude. 
And because they both scale with A 2 it would be impossible to separate the two by 
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performing grid refinement. Therefore it is critical 1o use a higher order accurate' 
numerical method when performing LES computations. 

Previous st udies [31-33] have applied LES codes to high Reynolds number jet 
flows. In each of the cited cases, agreement with experimental data has been poor. 
A common problem with these studies is the overprediction of the length of the jet's 
potential core. Each solution underpredicts the large scale structures and resultant 
turbulent mixing. In all three studies the computational domain began at the noz- 
zle exit and a “top-hat inflow velocity profile was specified. The nozzle geometry 
and influence of the jet lip was not modeled. Mankbadi et. al. [31] also artificially 
pert urbed the inflow conditions to in an attempt to excite the jet and increase mixing. 

1.4 Discussion of the Present Work 

The objective of the present work is to accurately predict a high Reynolds number 
turbulent nozzle flow field using computational methods and gain an improved under- 
standing of the jet's behavior. A completely new analysis code has been developed 
for this purpose. To meet the goal of accuracy, the sources of error present in any 
CFD analysis must be minimized. These errors, discretization error and modeling 
error are thoroughly investigated. 

Discretization error is addressed through the examination of an existing high-order 
numerical scheme and the development of a new scheme. Care is taken to consider 
both temporal and spatial accuracy conjointly. Both schemes are then examined 
analytically through their truncation error and experimentally through numerical 
experiments to ascertain their performance in solving representative problems. Grid 
refinement is also examined to show how increased resolution affects the solution. 


N AS A/TM — 200 1-210716 


8 



The largest contribution to the modeling error, modeling of the turbulent struc- 
tures, is examined fully. The Large- Eddy simulation technique was selected as the 
most promising method to use. 1 he LES equations are derived in detail and tin 
resulting sub-grid terms are examined. Models for the dominate sub-grid terms are 
implemented in the flow solver. Boundary conditions and other modeling issues along 
with correct code implementation are tested using simple validation cases and a rep- 
resentative jet problem. 

Using the knowledge gained from testing the codes on simpler problems, a tur- 
bulent compressible round jet is simulated. The jet has a Reynolds number of 1.2 
million and an exit. Mach number of 1.4. A complete description of the jet flowfield 
including turbulence information is presented. The time averaged flowfield is com- 
pared to experimental data to ascertain the accuracy of the simulation. Instantaneous 
flowfield data and turbulent statistics lend insight into the complex behavior of the 
jet. Correlation of velocity signals in the jet s mixing layer help quantify the large 
scale structures present. 
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Mixing Layer 


Figure 1.1: Jet schematic 


NASA/TM — 200 1-210716 


10 



CHAPTER 2 


GOVERNING EQUATIONS 

The equations governing the flow of the unsteady compressible jet are presented 
below. The fluid is assumed to be a continuum, which implies that the smallest scales 
of interest are much larger t han the scales of molecular mot ion. 

2.1 The Navier-Stokes Equations 

The most general set of equations considered are the three-dimensional Navier- 
Stokes equations. As presented below, they express the conservation of mass, mo- 
mentum, and energy for an unsteady compressible fluid in tensor form using cartesian 
coordinates (x,y,z). 

The continuity equation expresses the conservation of mass. 

9p + dpu 1 = Q (2.1) 

dt dxi 

The momentum equation expresses Newton's Second Law. It relates the time rate 

of change of momentum to the forces applied. 

OpUj dpUjXtj 0p_ _ 

dt dxj dxi dxj 
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The stress tensor is defined as 


Vij = - + 2//.S',-; 


and the strain rate tensor is 


' S - = 9 


1 / dit j da, 

2 V Ox, t)X; 


Sutherland's Law is used to model the viscosity 


f,7'> 


' c 2 + t 7 

The constants for // expressed in English units ( ) are C\ = 2.27 x L0~ 8 and 
C 2 = 198.6. 

Conservation of energy expresses the first law of thermodynamics. It relates the 
time rate of change of energy to the amount of heat added and the work done. 

dptt Opujet duiP _ dujajj _ % 

Ot dr i dx, dr, dr, 

The total energy of the fluid. e< is defined from the internal energy, t = c v T 


€t — e + 9 u i u ( 


The heat flux is represented by 


q, = -k- 


The system of equations is closed using the equation of state for a perfect gas. 


p = pRT 
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2.2 Axisymmetric Formulation 


A solution to the three-dimensional Navier-Stokes equations requires a large amount 
of computer memory, storage and C PU time. Taking advantage of the svmmetiv of 
a problem is a common technique used in ( FD to reduce the computer resources 
required for an analysis. Since the geometry of the round nozzle is symmetric about 
its centerline, it may be possible to model the problem using the axisymmetric form 
of the Navier-Stokes equations. In order for this assumption to be valid the flow field 
must also be symmetric about, the centerline. 

The axisymmetric form of the Navier-Stokes equations assumes that there aie 
no gradients or velocity components in the circumferential direction of a cylindrical 
coordinate system. The symbols r. and 9 represent the axial, radial, and circum- 
ferential directions. The corresponding velocities are represented by w, r, and u\ The 
resulting equations are presented below. 

Continuity 

ap + d £ u + [3pr r =0 (210) 

3t dx r 3r 


Axial momentum 


dpu 

dt 


+ 


dpu 2 1 dpruv 

dx r dr 



da, T 

17 + 


d<T xr 

dr 


+ <*xr 



(2.11) 


Radial momentum 


dpv 

~df 


+ 


dpuv 1 dprv 2 

dx r dr 



d(7 xr 

~d7 


dc T rr 

"I - X d - ® rr (XQ9 

dr 



2 d 

7 dr 
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The components of the stress tensor are 


<7r.r = f 
<Trr ~ f 
<7xr = ft 


°oo = ft 


j 4 du 

2 dv 

V3 d7- 

~ 3 lh- 

f A dv 

2 du 

\ 3 dr 

~ tlh- 


' du dv 
. dr dx 

2 / da dv \ 
2 V d.r dr J 


+ 


1 r 

3 7 


Energy 


dpt-t 

dt 


+ 


d/>t t u 

d.r 


1 dpe t rv dpu J dprv 8 

+ — H = T— {U(7 rx + V(J rr 

r dr d.r r dr d.r 

+ — (UC T xr + V(7 rr - q r ) + UVxr + V(7 rr 
dr 


rJi — 

3 r 


d ri v 1 
dr \ 3 ^ r 


d (2 uv 


The corresponding heat fluxes are represented by 


qs = 



q r 



2.3 Flux Vector Forms 


For OFD, it is helpful to recast the equations in flux vector form. This 
the equations the same form for easy discretization and solution. 

2.3.1 Three-dimensional 

The flux vector form of the three-dimensional X'avier-Stokes equations is 

OQ dE OF dC DE V dF v dG v 

dt dx Oy dz dx dy dz 


(2.13a) 

(2.13h) 

(2.13c) 

(2.13d) 

(2.14) 

(2.1oa) 

(2.15b) 

puts all 
(2.16) 
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where 


Q = 


p 

pa 

P v 

pw 

P*'t. 



pu 


0 


pu 2 + p 


ft XX 

E = 

pur 

E v = 

ft xy 


puiv 


<?r~. 


iptt + p) »_ 


ua xx + V(7 xy + wa xz — q x _ 


pv 


0 


pvu 


ft yx 

F = 

pv 2 + P 

Fr = 

ftyy 


pvir 


<T yz 


{ptt + p) 


+ I'ftyy H - l^ftyz m 


pw 


0 


pint 


ftzx 

G = 

pun ? 

(*V — 



piv 2 + p 


G zz 


{pet + p) u\ 


u(T zx + v(T zy + xr<T zz — q z _ 


2.3.2 Axisymmetric 


The axisymmetric equations are similar to the three-dimensional form with the 
addition of a source term. When the source term Ft is omitted, the equations simplify 
to the two-dimensional/planar form. 


f + ^ + SF + H= ^ + OK +Hi 

dt dx dy ox ay 


( 2 . 17 ) 


where 


P 


Q = 


pu 

pv 


ipet] 
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0 


r 

pil 


0 

l>a~ + p 

E, = 

°r.r 


puv 


(T. rr 

.(/* 

l + p) »_ 


U(T rr + t'(T rr - (Jj._ 

~ 

pv 


0 


pru 

F, = 

<*rr 

pi 

+ p 


(Trr 

.(/H 

t+v) v - 


tier,..,. + rcr,.,. - 



pr 


0 

,i = i 

V 

puv 

pv 1 

H t ,= 

a rr - (//*)_ 

a rr CTOS IfF §r£ {nL) 


+ p) 


ua x , + va rr - q,. - \p~- r— (%p£) - r§ (§ pf) m 
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CHAPTER 3 


NUMERICAL METHODS 

The governing equations are solved using two explicit finite difference methods, 
a MacCormack type predictor-corrector method and a Runge-Kutta method. Both 
types of schemes are in wide use today for fluid dynamic and acoustic analyses of 
nozzle flows. Each scheme has an associated set of strengths and weaknesses which 
will be discussed later in this document. 

The numerical schemes will be presented in terms of a one-dimensional model 
equation 

^ + ^ = 0 ( 3 . 1 ) 

dt dx 

where f = f(q). Extending the schemes to the Navier-Stokes equations (2.16 k 2.17) 
is a straight-forward matter. 

3.1 Predictor-Corrector Schemes 

MacCormack developed a two-step explicit finite difference method [34] which is a 
variant of the Lax-Wendroff scheme [35]. This method is very easy to understand and 
implement and hence has become very popular in the CFD community. The scheme 
is very robust and requires only two storage locations for each dependent variable 
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(2-N storage). As a consequence several variants of the MacCormack technique have 
been developed. The original technique and a higher order variant will be discussed. 

3.1.1 MacCormack 2-2 Scheme 


Mac( 'orniack s original scheme is second order accurate in both time and space. 
That is to say the truncation error of the scheme is proportional to both the time 
step and spatial step to the second power. 

The first stage, or predictor step, computes the solution at an intermediate time 
based on the solution at the previous time step. It uses a one-sided forward difference 
for the spatial derivative. 


9; = ( U 


At 

a7 


(/;+, - /;) 


(*■•*>) 


The second stage, or corrector step, computes the solution at the end of the time 
step based on the solution at the intermediate time. It uses a one-sided backward 
difference for the spatial derivative. 

The leading truncation error terms for MacCormack’s scheme are 

(A/) 2 d 3 q (A.r) 2 d 3 f At(Axf d 4 f 

6 dP 6 dx 3 24 dr 4 l ' " j 

where / has been linearized by / — Aq. The error associated with the temporal and 
spatial terms are dispersive in nature, while the error associated with the cross term 
is dissipative. This dissipative term is scaled by the time step. 




(3.2b) 


3.1.2 Gottlieb- Turkel 2-4 Scheme 


Gottlieb and Turkel [9] modified MacCormack's scheme to be fourth-order accu- 
rate in space while retaining the second-order time accuracy. They simply modified 
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the difference stencil of the spatial derivative to achieve higher accuracy. 


1 At 

6 A7 


(-/r +2 +s/r + ,-7/r) 


(3.4a) 


</f +1 = 3 


'll i * 

<li + ( U 


1 Ai t-r 

6~K~t ' ’ 


s/V. + !U) 


(3.4h) 


The leading terms in the truncation error for the (lottlieb-Turkel scheme are 

{At) 2 fPq (Ax) 4 ffq At (Ax) 2 c) 4 f , 

— + 18 9x 4 

In addition to retaining the second-order time accuracy, it is important to note that 
the cross term is dissipative and is scaled by At (Ax) . Bayliss et. al. [36] extended 
this scheme to the solution of the Navier-Stokes equations. 


3.2 Runge-Kutta Schemes 

As seen with the Gottlieb-Turkel scheme, all higher order variants of MacC’or- 
mack’s method maintain second-order time accuracy. It will be shown later that one 
cannot separate the temporal and spatial accuracy. They are clearly equated thiough 
equation (3.1). 

A true fourth-order accurate scheme in both time and space has been developed. 
A fourth-order central difference for the spatial derivative is combined with a fourth- 

order Runge-Kutta time stepping scheme. 

Runge-Kutta schemes are a popular family of numerical schemes with higher order 
temporal accuracy. These multi-stage schemes can be formulated for any order of 
accuracy. The number of stages in the scheme is equal to or greatei than the desiied 
order of accuracy. Two fourth-order Runge-Kutta schemes are investigated. 
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3.2-1 Standard Scheme 


The standard four-stage fourth-order scheme as given by .Jameson [37] is 


(ft) 

= </" 


</l 

= <r 

-l-A7D(^o) 

<72 

= q n 

-1a/D( 9i ) 

<73 

= <7 ’’ 

-1a7D(^) 

+1 

= q n 

— A/D (c/; 5 ) 


(3.6) 


The operator D is the spatial finite difference operator. For equation 3.1, D (q) would 
be a fourth-order finite difference stencil for — %}-. A central difference stencil for ~ 

dx I 

is used here. 


_ ~.fi+ 2 + 8./ i + l — 8/,_ 1 + .1 1-2 

dl- ~ ’ 12(A.r) 


(3.7) 


This scheme requires two storage locations for each dependent variable (2-N storage). 
The leading terms in the truncation error for this scheme are 


( At) 4 d 5 q (Ax) 4 _ A/(A.r)-' ^8 6 f 


120 dt° 30 da 

3.2.2 Low Dispersion Scheme 


30 8:r 


(3.8) 


Several researchers [38-41] have developed alternative Runge-Kutta schemes that 
have a lower dispersion error than the standard scheme leading to greater stability 
and accuracy. To accomplish this, additional stages are required. The additional 
stages provide a means to impose the additional constraints necessary to minimize 
the error. All of the schemes are based on a general M-stage 2-N storage formulation 
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given by 


( iq m = 1 + A/D ( ({m - 1 ) (o.^a) 

</m = ( lm- 1 $md(fm (3..)b) 

for m = 1 ...A/, and where q u = r/ n and r/A/ = <y n+1 - The coefficient o-j is typically 
set. to zero for the algorithm to be self-starting. Again, the operator D is the spatial 
finite difference operator. 

Carpenter and Kennedy's five-stage fourth-order scheme [38] was chosen for its 
fourth-order accuracy, low number of stages, and ease of programming. The coeffi- 
cients for the scheme are given in table 3.1 

stage (m) oc m 

1 0.00000000000 0.1496590219993 

2 -0.41789047450 0.3792103129999 

3 -1.192151694643 0.8229550293869 

4 -1.697784692471 0.6994504559488 

5 -1.514183444257 0.1530572479681 

Table 3.1: Coefficients for fourth-order low-dispersion Runge-Kutta scheme 


The leading terms in the truncation error for this scheme are 

(At ) 4 cfq (Ax )*d*q At(Ax)\(ff ,, im 

"3oTaF + _ M _ ac so A rh 6 

The dispersive error term due to the time step (the first term) is two and one-half 
times smaller than the standard fourth-order scheme’s error (3.8). 
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3.3 Fourier Analysis of Numerical Schemes 


In the field of high order numerical methods for DNS and EES. the proper res- 
ohttion of the waves present in the flowfield is rritiral. A popular met hod to assess 
the ability of a numerical scheme to resolve waves is the Fourier analysis of the semi- 
discretized equation (spatial discretization only) [42]. 

The semi-discretization of equation 3,1 is 


dq " n i ■ 


(3.i i ; 


where q„ is the discrete solution. We choose a. sinusoidal trial solution with wave 
number u\ which has the form 


<?„(-'. /) = (3.12) 

w e substitute 3.12 into 3.11 and solve for v(u,\t). The resulting solution of the 
equation is 

q n (u,\t) = v{uj. o)g Re [ £, ^ )t ] e *"( J ’>-^d (3.13a) 

with 

c‘ = -Jlm[D(u?)#] (3.13b) 

where D(us) are the eigenvalues of D. The exact solution to equation 3.1 is 

q{w,t) = v{u;.0)e MT - ct) (3.14) 

Bv comparing the numerical and exact solutions we can see how the numerical scheme 
affects the propagation of the wave. The exact solution (equation (3.14)) shows that 
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amplitude of the wave should not change. However, in the numerical solution the 
amplitude is a function of time (equation (3.13a)). Therefore. 


if Re [/)(u?)] = 0. then the scheme is conservative 

if Re [£)(u;)] < 0. then the scheme is dissipative 

if Re [l)(u;)] > 0. then the scheme is unstable 


Similarly, comparing the numerical wave speed c* to the exact wave speed c we find 


if — = 1, then the scheme has no phase error 
c 

if — ^ L then the scheme introduces a phase error 
c 

To achieve a conservative scheme (Re [D(u;)] = 0) the matrix D must be anti- 
symmetric. Central difference schemes, such as those used here, satisfy this criteria. 
Upwind schemes, which are dissipative, do not have an anti-symmetric mattix and 
were not considered suitable for this study. The wave speed relative to the exact 
speed for each scheme considered here is 


c* 

c 

c* 

c 


sin(u?Ax) 
uj Ax 

-^sin(^A.r)co.s(^’Ax) + ±sin(u:Ax) 
ojAx 


. for second-order spatial accuracy 

(3.15a) 

. for fourth-order spatial accuracy 

(3.15b) 


Figure 3.1 compares the error in the wave speeds for second- and lourth-order 
accurate central difference operators, D. The error is a function of the wave numbei. 
u;. and grid resolution. Ax. From this Fourier analysis perspective the numerical 
scheme can be seen to act as a spectral filter of the exact solution. This enoi is 
sometimes plotted as numerical wave number versus exact wave number (/Ar versus 
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a.’A.r. figure 3.2). It is dear that higher order schemes can resolve higher wave number 
waves with less grid resolution. But , neither scheme can adequately resolve very high 
wave numbers, even with very fine grid resolution. 

If one determines an acceptable error in wave speed for their calculation, the 
maximum value ol u,’A.r that provides t his error can be related to the number of grid 
points per wavelength necessary to accurately resolve a wave [43]. 


N = 


h r 


(3.16) 


The wave will be properly resolved if it is captured by at least N points. This "points 
per wavelength" has become a popular measure of a scheme's "goodness". But, this 
measure does not consider the error due to time discret ization and the computational 
cost of the scheme. 


3.4 Filtering 


Both numerical methods introduce a dispersive error into the solution. This er- 
ror manifests itself in two ways. First, as shown above, waves that are that are not 
resolved with enough grid points have their speeds altered. Secondly, new high wave 
number waves are introduced into the solution. If unchecked, these errors can grow 7 
and cause the analysis to become unstable. To damp out these waves, an artificial 
dissipation term is typically added to the equations [44-47]. All these ^classical" 
techniques involve the use of a user specified constant that scales the effect of the dis- 
sipation. The value of this constant is somewhat arbitrary and varies widely between 
cases. 

In this study a solution filtering technique is used. This technique can be regarded 
in two ways. In the classical sense, solution filtering simply adds a dissipative term to 
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the equations to damp the unwanted waves. From a Fourier analysis perspective, this 
technique filters the high wave number components out of the solution. A properly 
selected filter removes both the numerically introduced waves and the poorly resolved 
waves, leaving only the portions of the flowfield that are accurately modeled. Sev- 
eral researchers have developed solution filters for ('FD including [10.42. 48.49]. The 
explicit filters of Kennedy and Carpenter [48] were chosen for their ease of implemen- 
tation and clear relation to the unfiltered solution. The filtered solution is simply t he 
unfiltered solution plus a dissipative term. 


</< 


</, — o/j( A.r) 


2 71 


d 2n qj 

d.r 2n 


(3.17) 


where oo — — ^ — and n — — 1.2,. .. . Kennedy and Carpenter developed a familv 
of filters with corresponding boundary stencils for n = 1...7. These filters are 
implemented in the analysis code. Bv choosing the order of the filter. 2ft. to be larger 
than the order of accuracy of the numerical scheme, one can insure that the filter 
does not overly influence the numerical solution. The response of the filters is shown 
in figure 8.3. As the order of the filter is increases the cut-off wave number of the 
filter is increased and a greater portion of the domain remains unmodified. 

Examples of the effect of solution filtering are shown in figures 3.4 and 3.5. A 
sine wave on the domain 0 < .r < it was filtered 100.000 times (a typical number of 
iterations for a jet calculation). To illustrate how increased grid resolution reduces the 
effect of the filter, sine waves made up of 5, 15, 25, 35, and 45 points (u;A.r = 1.257, 
0.4189, 0.2513, 0.1795, and 0.1396) were filtered using fourth-, sixth-, and eighth- 
order filters. Figure 3.4 shows that increasing the resolution of the wave reduces the 
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effect of the filter. Poorly resolved waves can be completely eliminated (the fourth- 
order filter completely removes the 5 point resolved wave) Figure 3.5 shows how t he 
resolved wave number increases with the order of the filter. A wave resolved by five 
points was filtered with fourth-, sixth-, eighth-, tenth-, twelfth-, and fourteenth-order 
filters. The wave (u,’A.r — 1.257) is completely removed by the fourth order filter, but 
damped less than 10 percent by the tenth-order filter, and maintained exactly bv the 
fourteenth-order filter. 
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Figure 3.1: 


Error in wave speed 
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to Ax 


Figure 3.2: Error in wave number 
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X 


(a) fourth-order filter 



x 


(b) sixth-order filter 


Figure 3.4: Effect of grid resolution on filter response 

continued 
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CHAPTER 4 


1-D ERROR ANALYSIS 


The behavior of the numerical schemes presented in Chapter 3 is examined for 
three one-dimensional model problems. Truncation error, efficiency, and the conse- 
quence of disparate temporal and spatial accuracy are discussed. The schemes are 
used to solve the one-dimensional inviscid convection of a gaussian pulse. 

The accuracy of these schemes is commonly expressed separately in terms of spa- 
tial and temporal accuracy. A great deal of effort has focused on increasing the spatial 
accuracy of a numerical scheme without regard to the temporal accuracy. I his chap- 
ter examines the behavior of two commonly used schemes in terms of truncation error 
and computational cost on a model equation. Particular attention is paid to the trun- 
cation error of the schemes and how the spatial and temporal errors affect the overall 
order of accuracy of the scheme. 

Three one-dimensional problems were used to test the accuracy and efficiency of 
the schemes. All three are based on the model equation (3.1). 
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4.1 Problem Descriptions 

4.1.1 Problem 1 


The first problem is the linear convection of a gaussian pulse, where q = / = u. 

(4.T 


<ju chi 

¥ + 57 = 0 


The initial condition is given by 


«(*,<)) = u 0 (;r) = i f - ,n(2 >(§) i 


(4.2) 


The domain is —20 < x < 450 and the solution is run for 0 < t < 400. An exact 
solution to this problem exists and is given by 


u fX (x, t) = u 0 (x - t) 


(4.3) 


4.1.2 Problem 2 


The second problem is nonlinear where q = u and / 

d u c) ( u 2 


L y- 


m + rAj ) = 0 


4.4) 


The initial conditions and simulation time are set so that a smooth final solution is 
obtained. The initial condition is 


u(x, 0) = u 0 (x) = |e in(2) (io) 


(4.5) 


The domain is —50 < x < 50 and the solution is run for 0 < t < 100. A numerical 
approximation to the exact solution is obtained using the method of characteristics. 
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4.1.3 Problem 3 


The third problem is also nonlinear (equation (4.4)). but the initial conditions and 
time of the simulation are set so that the pulse is allowed to coalesce into a shock. 

The domain is —50 < x < 50 and t he solution is run for 0 < t < 200. Like problem 
2, a numerical approximation t.o the exact solution is obtained using the method of 
characteristics. The location of the shock is then fitted in the exact solution using 
Whitham’s area rule [50]. 

4.2 Results 


All the one-dimensional calculations were run on an Apple Macintosh Powerbook 
Cl 3 computer with a 266MHz PowerPC G3 processor. 

The standard scheme was used for all the Runge-Ivutta results in this section. 
The error of the numerical scheme was measured by the I 2 norm, which is computed 
as follows 


lo = 


rin V ' / \2 

/ M l i ~ ) 


(4- 


l=\ 


For each problem, both numerical schemes were run to determine the maximum 
stable time step. Then, each scheme was run at a number of different values of as 
the spatial step was halved. 


4.2.1 Problem 1 

A sample solution to the linear problem for both numerical schemes is presented in 
figure 4.1. The Runge-Kutta scheme shows better agreement with the exact solution 
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t han the Gottlieb-Turkel scheme. Error from the Gottlieb-Turkel scheme is plotted 
versus spatial .step and time step (figure 4.2). The terminal slope of the line, p 
indicates the order of accuracy of the scheme and is listed in the plot legend. The 
slopes in figure 4.2(a) show that the error is of second-order accuracy for t he varying 
spatial step. Only where the time step is much smaller than the spatial step does the 
scheme behave as a fourth-order scheme. Also, the error is dependent on the time 
step chosen. This is due to the third term in the truncation error (equation (2.5)). 
Clearly, this dissipative error is significant. Figure 4.2(b), shows the error behavior 
with the time step. Again, the scheme is only second-order accurate except when the 
time step is very small compared to the spatial step. In order to obtain fourth-order 
behavior from the Gottlieb-Turkel scheme, the error due to the time step must be 
reduced at the same rate as the error due to the spatial step. In order to accomplish 
this, the time step must be reduced by a factor of j as the spatial step is reduced by 
y. Figure 4.3 verifies that when adjusting the time step in this manner fourth-order 
accuracy is obtained. 

Results for Runge-Kutta scheme, on the linear problem (figure 4.4), indicate that 
the scheme is truly a fourth-order accurate scheme. In addition, except for the highest 
value of ^ where the time error term dominates, the error is independent of time 
step. 

The efficiency of a numerical scheme can be seen by comparing the l 2 error against 
the time required to obtain that error level [51]. Both schemes are compared at their 
maximum stable time step. Figure 4.5 shows that the Runge-Kutta scheme is superior 
to the Gottlieb-Turkel scheme. Runge-Kutta achieves an equivalent error level in 
about an order of magnitude less time than the Gottlieb-Turkel scheme. This result 
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is due to the lower truncation error and larger allowable time step of the Runge- 
Kutta scheme. Reducing the time step of the Gottlieb-Turkel scheme to maintain 
fourth-order accuracy increases the accuracy but also increases (he computational 
cost. 

4.2.2 Problem 2 

Sample solutions to the smooth nonlinear problem for both numeiical schemes Hie 
presented in figure 4.6. Because a short simulation time was necessa.iv to maintain a 
smooth solution, the differences between the exact solution and numerical solutions 
are not discernible on the plot. Error data for the smooth nonlinear problem are shown 
in figures 4.7- 4.10. The effects of round-off error, the reduction in accuracy/slope at 
small A/ and A.r. make it difficult to report a terminal slope, instead the maximum 
slope was used. The trends for the nonlinear analysis are similar to the linear case. 
However, the benefit for using the Runge-Kutta scheme is somewhat less. But, at 
large time steps, for maximum efficiency, the Runge-Kutta. scheme is superior. 

4.2.3 Problem 3 

Solutions to the nonlinear problem with a shock (figure 4.11) show that both 
schemes produce high frequency oscillations near the discontinuity. This problem does 
not provide useful error information and illustrates the limitations of finite difference 
schemes. The presence of the discontinuity reduces the accuracy of all schemes to 
first-order (figures 4.12 & 4.13). Results from the Gottlieb-Turkel scheme indicate 
that the third term in equation (3.5) provides a damping effect, on the numerical 
oscillations that occur near the shock. This damping allowed the scheme to obtain a 
converged answer. Figure 4.12 shows that larger time steps yielded lower errors for 
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this problem. The dissipative term in the truncation error is scaled with the time step 
and larger time steps resulted in greater dissipation of the non-physical oscillations. 

To obtain useful solutions to this problem, the Runge-Kutta scheme required 
some form of artificial dissipation. The solution filtering technique of Kennedy and 
Carpenter [48] was used. Filtering adds to the solution an additional dissipative term 
that removes high frequency numerical oscillations while leaving the low frequency 
physical oscillations untouched. 

The order of the filter. 2 n in equation (3.17). determines the magnitude of the 
dissipation added and the range of frequencies which are damped. If the order of the 
filter is larger than the order of the truncation error of the numerical scheme, then 
the filter should have a negligible effect on the solution error. Increasing the order of 
the filter increases the range of frequencies which are untouched by the filter. Figure 
4.13 shows that the filters remove the non-physical oscillations from the solution and 
greatly reduce the error of the scheme. The insensitivity of the error to the order of the 
filter indicates that the truncation error of the scheme is greater than the dissipative 
term of the filter and that the numerical oscillations have frequencies higher than 
frequency range of the highest order filter. 
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(a) Spatial error 
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(b) Temporal error 


Figure 4.2: Error of Gottlieb-Turkel scheme for the linear problem 
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Ax 


Figure 4.3: Improved spatial error of Gottlieb-Turkel scheme for the linear problem 
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(a) Spatial error 
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(b) Temporal error 


Figure 4.4: Error of Runge-Kutta scheme for the linear problem 
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Figure 4.5: Efficiency of the numerical schemes for the linear convection problem 



Figure 4.6: Solution to smooth nonlinear convection problem 
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(a) Spatial error 
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(b) Temporal error 


Figure 4.7: Error of Gottlieb-Turkel scheme for the smooth nonlinear problem 
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Figure 4.8: Improved spatial error of Gottlieb-Turkel scheme for the smooth nonlinear 
problem 
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Figure 4.10: Efficiency of the numerical schemes for the linear convection problem 
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(a) Complete wave 



35 36 37 38 39 40 41 


x 


(b) Detail of shock 


Figure 4.11: Solution to nonlinear convection problem with shock 
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(a) Spatial error 



(b) Temporal error 


Figure 4.12: Error of Gottlieb-Turkel scheme for the nonlinear problem with shock 
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(a) Spatial error 





(b) Temporal error 


Figure 4.13: Error of Runge-Kutta scheme for the nonlinear problem with shock 
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CHAPTER 5 


LARGE-EDDY SIMULATION 

In this chapter the equations for Large-Eddv Simulation are presented. The mod- 
eling of the unresolved (sub-grid) terms and implementation into the flow solver is 
discussed. 

The basis of LES is the separation of the large and small scale turbulent fluctu- 
ations. To separate the large (resolved) and small (unresolved) scales the equations 
are filtered. A spatial filter G with a filter width A is used. 

/= - amw (b.n 

The overbar represents the resolved, filtered, or large-scale portion of the function. 
Commonly used filter functions are a box filter, a Gaussian filter, or a spectral cutoff 
filter [30]. Typically, it is not necessary to know the exact form of the filter function 
G, but only that it exists. In practice, the solution is not explicitly filtered. It is 
assumed that the numerically computed solution is a filtered representation of the 
exact solution. This assumption is justified based on the Fourier analysis results of 
sections 3.3 and 3.4. There, the numerical schemes and solution filtering were shown 
to behave as a spectral cutoff filter of the exact solution. 
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Several constraints are imposed on the filter function. 


1) G(-0 = <7(£) 

2 ) fl G ^ )d z = 1 

3) <7(0 ->0 a« |^| — ► oo 

4) 6'(0 is small outside { — j- j) 


These constraints are necessary to insure that the filter function will commute with 
the derivative. 


<!L = d L 

dx dx 


(5.2) 


Favre (density) weighting is used in the filtering process. This allows for convenient 
recovery of terms corresponding to the unfiltered equations. 

/ = 2 
P 

5,1 Filtered Equations 


The filtering process is applied to the continuity, momentum, arid energy equa- 
tions. Details of this process can be found in Appendix A. The resulting equations 
are comprised of resolved and unresolved terms. The resolved terms in the filtered 
equations directly correspond, in form, to the unfiltered equations. The additional 
unresolved or sub-grid terms are modeled as source terms to the equations. 
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5.1.1 Continuity Equation 

The filtered continuity equation contains only resolved terms. Thus, it corresponds 
direct lv to the unfiltered equation (2.1). 


dp = 0 

f)t ().r, 


(5.4) 


5.1.2 Momentum Equation 


The momentum equation contains two unresolved terms (underbraced) that must 
be modeled. 


dpuj , dpujUj , dp da t j drij d ~ ] 

~df + ~07~ + d7~ rlr, dx 3 + dxj [ U ' J ‘ 


Three stress tensors occur in (5.5). The filtered stress tensor 




a l3 = -§ pSijSkk + 2 fiSij 


(5.6) 


The resolved stress tensor 


a tJ = — | p-SijSkk + 2pSi j 


where p = //(T). and 


1 / duj 8u; 
S; . = -1-^+- 


(5.8) 


‘~ tJ 2 V dxi ' dxj 

The third is the unresolved or sub-grid scale stress tensor, r tJ . It is this term that 
provides the effect of the sub-grid scale turbulence on the larger resolved scales. 
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5.1.3 Energy Equation 


The Favre filtered energy equation contains five additional unresolved 
derbraced) that must be modeled. 


dpT t dpiiitt duip dujdij 


dt 


rh 


Ox, 


Ox, 


()<), o_ 

().v, dx, 


i J + f — 


dpDi 

dxi 


where the filtered total energy is 


d ( do a _ ddn\ 0 

+ 77— I «> — - u,-r- I - - 7 — (q t 


dxi \ J dx l J dr, ) dx, 


e 1 : e + 2 UkUk 


and the filtered and resolved heat flux vectors are 


<h = —k 


OT 


dx, 


t dT 

= -k — 

ox. 


where k = k(T). The sub-grid scale heat flux is given by 


Qi = p (uiT - 


The sub-grid scale turbulent dissipation rate is 


Oil , „ diij 


arid the sub-grid scale turbulent diffusion is 


pD{ = | (wu- - UiUkUk) 


terms (un- 

q,) (5.9) 

(5.10) 

(5.11a) 

(5.11b) 

(5.12) 

(5.13) 

(5.14) 
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The equation of state becomes 


p = pR'f (5.15) 

Since the total energy, equation (5.10). contains the unknown Uk»k- pressure must be 
obtained using the following expression 

p = (7 - 1 ) (p<~t - - |m.) (b.J 6 ) 

5.2 Sub-Grid Scale Modeling 


The unresolved or sub-grid scale contributions to equations (5.9) and ( 0 . 0 ) must 
be modeled. The methodology used here is based on the mcompiessible sub-giid 
scale model of Smagorinsky [24]. Additional terms to account for compressibility 
were added based on the work of Moin et. al. [2b] and \ reman et. al. [•)-]• 


5.2.1 Momentum Equation 


The Smagorinsky model is a very popular and widely used model foi the sub 
grid scale stress tensor. It. is an eddy viscosity model where the sub-grid scale stress 
tensor is modeled as an eddy viscosity multiplying the resolved stress tensor. It was 
developed for incompressible flows and has been frequently extended to compressible 
flows. The compressible form of the model, given by Moin et. al., was used here. 

m = §f / pA' 2 |A| 2 d u - 20A i |5| {Sij - ±S kk 6 tJ ) (5.17) 


The coefficients C and Cj are user defined constants. The second unresolved term in 
the equation results from the nonlinearity of the viscous stresses and is neglected. 
The final form of the modeled momentum equation is 


6pu, dpuiUj dp _ d 

dt + dxj dxi dxj 


ft 


(5.18) 
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The eddy viscosity is 


//, = CpA' 2 \S\ ( 5 . 10 ) 

and <I> is 

<D = §C/pA 2 |5j 2 (5.20) 


5.2.2 Energy Equation 


There are many different sub-grid models for the energy equation and no one 
set of models is as popular as the Smagorinsky model for the momentum equation. 
This is because there are many forms of the energy equation that result in different 
sub-grid scale terms. It is also because the majority of work in LES has been done 
for incompressible flows. The filtered energy equation ( 5 . 9 ) contains five unresolved 
terms. 


The first term, the sub-grid scale heat flux is modeled based on Moin's work [28] 


Qi 


^ dT 
Pr t Ox, 


(•5.21) 


The model for the second term comes from Vreman [52] and is given by 


£3 

2A 


Tkk 


(5.22) 


The coefficient f 3 is a user defined constant. The last three terms are typically 
considered to be much smaller than the other terms and are neglected. 

The final form of the modeled energy equation is 


Opt, dpu,7, duip 0 

dt dxi + dx{ dxi 


Uj(Tjj 4- I p + // f 


Pr\ 7 R 1 ST' 


Pr,J 7 -I Prdxi 


+ c (5.23) 
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5.3 Implementation 


The sub-grid scale terms in equations (5.18) and (5.23) and their correspond- 
ing models were implemented in both the two-dimensional/axisvmmetric and thiee- 
dimensional Navier-Stokes flow solvers. It is important to note that the modified 
definition of the pressure (equation (5.16)) must be used. 

Because the resolution of scales in a numerical scheme is directly related to the 
grid resolution, the filter width A was chosen to be a characteristic length of the 
computational grid. Since the grid is not uniform, this length varies widely over the 
grid. The filter width at each location was defined as the cubed root of the volume 
associated with each grid point. 

To maintain a laminar sub-layer in wall bounded regions, the effect of the sub-giid 
model must be diminished near the wall. To accomplish t his, a Van Driest damping 
funct ion [53] was used to scale the effect of the sub-grid terms. 

/v D =( (V2.1) 

where y + is the inner variable distance and the constant A + is set to 26. 

The constants for the sub-grid models were chosen based on previous studies. 
Erlebacher et. al. [54], Moin [28], and Vreman [52] derived the model constants based 
on DNS results. Erlebacher concluded that, C = 0.012, and Moin gave a range of 
values of, 0.008 < C < 0.014. Moin also provide a. range of values for C/, where 
0.0025 < Ci < 0.009. From these data the coefficients were chosen to be C - 0.012, 
and Ci = 0.00575. Vreman determined that C 3 = 0.6. 

Results from LES simulations are obtained in terms of filtered conservation vari- 
ables. These variables differ from the exact conservation variables by equation (5.1). 
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However, since the form of the filter function. G, is unknown, we cannot recover the 


exact conservation variables. We must, assume that the difference between the two 
forms is small. For convenience, when reporting LFS simulation results the overbar 
( ) and tilde ( ~ ) are dropped in the notation and must be assumed. 
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CHAPTER 6 


FLOW SOLVER 


The numerical methods outlined in Chapter 3 are implemented in two computer 
codes to solve the unsteady Navier-Stokes equations. Separate three-dimensional and 
two-dimensional/axisymmetric codes were created. I he codes were written in the 
FORTRAN 77 computer language and should run on any computer platform with a 
FORTRAN compiler. 

While the purpose of the codes is to simulate compressible turbulent jets, they were 
written so that they can be easily adapted to a wide range of flows and geometries. 

Large-Eddy simulations can be performed by solving the filtered unsteady Navier- 
Stokes equations with the sub-grid scale model, or direct solution of the equations 
without approximation can be made. By neglecting the viscous terms in the Navier- 
Stokes equations the Euler equations can also be solved. 

6.1 Generalized Curvilinear Coordinates 

To allow easy solution to a wide range of geometries the Navier-Stokes equations 
are solved using generalized curvilinear coordinates [55]. This enables the computa- 
tional grids to be fitted to complex shapes and allows for grid stretching. 
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6.1.1 Coordinate Transformation 


Lasy and efficient implementation of the numerical scheme is done by working in 
a computational domain. I his domain is rectangular and consist s of equally spaced 
grid points. Tin- computational domain (£. //, () is mapped onto the physical domain 
(.r. y. ~) using the following transformation 


£ = £(•''• y— ) 

>1 = '/(•''•'/• = ) 

C ={(x,y.z) (6.1) 


Derivatives in carlesian coordinates are computed using the chain rule 


d_ _ d d d 

dx ~ + rh % + 

d _ c () d . d 

lfy-^ + th d^ + ^ 


d 

d~ 


d 


d 




(6.2) 


The terms f T . £ y , ij T , r/ v , //., C- Cy- and (\ are the grid metrics associated with 
the transformation. The grid metrics are computed numerically using the following 
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method. 


£,x — J{ih )~c Vc ~ 1 1 ) 

Sy = ■/ ( -l\ ~ c) 

6 = 

ih = •/(</<-€ - %~c) 

h.v = J(' v t~C — x c z 0 

n~. = J(xcyt-*tVc) 

Cr = Vv-i.) 

Cjr = J ( x n~£ 

C- = Ji^ejr, - x >)ye) 

The Jacobian of the transformation. J. is defined as 

J = ( 6 . 4 ) 

x i(yv z C - Vc z n) - -‘'riiyt-c - yc z s) + x <(fk z n - -In ~i ) 

The terms x^. x v , x^. y^. y n - (/<;. c rj . and are computed using fourth-order central 
differences in the interior of the grid and one sided differences on the boundaries. The 
calculation is easily done in computational space because the grid is evenly spaced in 
the »?. and ( directions. 

6.1.2 Chain Rule Form of the Governing Equations 

The chain rule form of the governing equations is 

^(.^E-E.n^E-K) + ^(E-E.) 

+ (F - F„) + .,»T {F - FJ + (F - F„) (6.5) 

4- (G - G,; + y — (G - G,.) -h (. jj- (G - G, ) = 0 
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This form of the equations is a weak conservation form. Most (TI) solvers use a 
strong conservation law form of the equations developed by Vinokur [56]. However. 
Hixon et. al. [57] recently showed that in practice, the chain rule form of the equa- 
tions is more accurate than the strong conservation form when the metric terms are 
computed numerically. 

6.2 Time Stepping 


The time step used to advance the solution is based on the relation 


At = CFL 



(6.6) 


where .4 is the local wave propagation speed and CFL is the Courant Friedrichs 
Lewey number [8], a constant which is based on stability considerations. The term 


A.r 


4 is a characteristic time and is computed using the inviscid CFL condition [58]. 

— : : — \ -l 


At 


CFL 


Ax 

“T 


U 


+ 


|u 


+ h a 


A.r Ay Ac 


1 


Ax 


A + + 


A y 2 A 


~2 


(6.7 


For time accuracy all points in the domain must be advanced at the same time step. 
Therefore. A tcFL is computed at every grid point and the minimum value is used to 
compute the time step. 


At = CFL - Aten 

6.3 Treatment of the Viscous Terms 


( 6 . 8 ) 


The computation of the derivatives within the viscous fluxes ( equations (2.4), 
(2.8), (2.13), & (2.15a) ) is done differently for both numerical methods. The predictor- 
corrector schemes require a fairly complex treatment of the viscous terms. Bayliss et. 
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al. [36] proposed the following method to maintain the order of accuracy of the Mac- 
Oormack type schemes. If the derivative within the viscous flux is being differenced 
in the same direction as the flux is differenced, then the derivative in the viscous flux 
is computed using a one-sided difference in the direction opposite of direction the 
flux is being differenced. If the derivative within the viscous flux is being diffVieiuf d 
in a direction other than the direction the flux is differenced, then the derivative in 
the viscous flux is computed using a central difference. Consequently, the derivatives 
in the viscous terms are computed twice (one-sided and central differenced) foi each 
step. The Runge-Kutta based schemes simply use fourth-order central differences for 
all viscous derivatives. 

6.4 Boundary Conditions 

Boundary conditions can be specified on any portion of any giid surface within the 
computational grid. This allows for a large range of complex shapes to be modeled 
within a single computational domain. Internal portions of the grid that aie enclosed 
by computational boundaries are specified as "holes”. The conservation variables 
within the hole regions are not updated during the solution process. Figure 6.1 
illustrates the use of internal boundaries and holes to model a nozzle wall. 

6.4.1 Inflow and Outflow Boundaries 

The formulation of these boundary conditions is based on the local one-dimensional 
propagation properties of the flow [46]. These properties are obtained from the eigen- 
values of the Jacobian matrix in the quasi-linear one-dimensional Euler equations. 
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The flow is assumed to he one-dimensional in the direction normal to the boundary. 
The resulting characteristic variables are (/,/?,. </,?/, + a. and </,/>, — a. when' it; is the 
unit normal vector to the boundary and a is the speed of sound. The variable v;n t 
occurs twice in two-dimensional flows and three times in three dimensional flows. 
These eigenvalues determine the propagation direction of the characterist ic variables 
at the boundary. 

Subsonic Inflow 

For a subsonic inflow boundary u, < a. This results in four (three in two- 
dimensions) positive characteristic variables and one negative variable. Therefore 
four variables propagate from the boundary into the domain and one variable prop- 
agates out of the boundary from the domain. To correctly mimic this behavior in 
the computation, so that the problem is mathematically well posed, four of the five 
conservation variables are specified at the boundary and one is determined from the 
interior. 

The flow from the inflow boundary is assumed to be normal to the boundary 
so that the transverse velocity components are set. to zero. The total pressure. p 0 . 
and total temperature. To are specified. The velocity normal to the boundary is 
extrapolated from the interior. For a constant £ surface (constant i ) with its normal 
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vector aligned in the x direction, the conservation variables are are 


obtained as follows. 


Po 


i 

7-1 


RTn 1 1 + - ll 


2 


't — 1 ( t + 1 j.k) 


(t+Fj.AO 


Q 


= o 


ij-h) 

Q*,,,, = 0 

Q r) U,j,k) 


RTo + 


7 - J 


(6.9) 


Supersonic Inflow 


For a supersonic inflow boundary it, > a. This results in five (four in two- 
dimensions) positive eigenvalues. Therefore all information propagates along the 
characteristics from the boundary to the interior. 

All five conservation variables are specified on the inflow boundary and are held 
fixed. For a constant £ boundary the conservation variables are 


Q'l{i ,j,k) P 

Q% Pc o^oo 

Poc W oc 

Qs { ,.j.k, = ^- 10) 

Subsonic Outflow 

For a subsonic outflow boundary it, < a. This results in four (three in two- 
dimensions) negative eigenvalues and one positive one. Here information from four 
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of the characteristics propagates from the interior to the boundary and information 
from one characteristic propagates from the boundary to the interior. 

The stat ic pressure. . is specified on the outflow boundary. Density and all three 
velocity components are extrapolated from the interior. Total energy is computed 
using the specified static pressure and the extrapolated density and velocities. 


t [ nnax ,j <k) 
tnui.v,j,k ) 

Q r 

tnia.T ,j .K) 

Q ^{tinux ,j ,k) 
imaij,/:) 


Q 1 ( I 771 « X — 1 , J . k ) 

Ql(n> lU'i-l.j.k) 

^(imax— 1 ) 

n Ql +QI + QI 

rt , ^ *{imax,j,k) * J (i max.j,k) ^ *(trnax,j. 

i 


k) 


7 - 1 


2Q h 

V 1 { t 771 ax,j 


(6.1 1 ; 


k) 


Supersonic Outflow 


For a supersonic outflow boundary u ( > a. This results in all five (four in two- 
dimensions) eigenvalues being negative. Information from all characteristics propa- 
gates from the interior out through the boundary. 

All fi ve conservation variables are extrapolated from the interior onto the bound- 
ary. 


1 ( tmax.j'k) 

** 1 { 1771 ax — 1 ,J ,k ) 

Q’2{i rnaxj.k) 


imax,j f k) 

- O'X 

* *^( trnax — 1 ,J ,k) 

Q^{trnax :J: k) 

iTTiaX — 1 ,J,k) 

Q^(imax,j,k) 

^(imar- 1 ,j,k) 
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Non-Reflecting Exit Zone 


Waves that propagate from the interior of the domain to the outflow boundaries 
can reflect off the boundary and back into the computational domain. These reflected 
waves are non-physical and may contaminate or distort the solution. These reflections 
occur because the extrapolation of interior information to the boundary is done along 
grid lines and not in the direction of the propagat ing waves. 

Expanding the computational domain so that the boundaries are far away from 
the area of interest reduce this problem, but adds significantly to the computational 
cost. Several authors have attempted to eliminate these reflections through the use 
of characteristic boundaries [59-61]. However, these efforts have met with limited 
success for two- and three-dimensional applications. 

A relatively simple and more effective method to reduce reflections is the use of 
"exit zones" [18,61,62]. In this method a simple outflow boundary condition, such 
as those outlined above are used. In the region adjacent to the outflow boundary, 
a combination of grid stretching and solution filtering is used to damp the outgoing 
and reflected waves. Rapid grid stretching increases A.r and therefore increases the 
associated dissipation in the truncation error of the numerical scheme (equations 
(3.5). (3.8) and (3.10)). Solution filtering (section 3.4) provides additional dissipation 
that damps a large range of wavelengths. 

6.4.2 Solid Surfaces 

Two options exist to model solid surfaces. Surfaces in viscous calculations are 
specified with a no-slip wall boundary condition. Surfaces in inviscid calculations 
and planes of symmetry are specified using a slip wall boundary condition. 
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No-Slip Wall 


The no-slip wall boundary imposes no relative motion between the fluid and the 
solid surface. Enforcing no-slip at the wall is done by setting all velocity components 
on the boundary to zero. The density and total energy are determined using two 
additional constraints. The first constraint stems from boundary layer theory. The 
pressure gradient normal to the wall is set to zero. The second constraint assumes no 
heat transfer through the wall, and consequently the temperature gradient normal to 
the wall is set to zero. For an ?/ constant wall (constant j ) . the conservation variables 
are written as 




Qi 










= 0 


= 0 




Ql + Q\ + Q\ 


1 (*.j+ 1 ,A) 

This formulation assumes that the grid lines are normal to the wall. 


( 6 . 13 ) 


Slip Wall 


An inviscid wall surface and a symmetry plane share the same constraints and 
hence the same boundary condition. The velocity vector of the flow along the surface 
must be tangent to the surface. To accomplish this, the velocity vector one point 
off the wall is decomposed into components that are normal and tangential to the 
surface. The normal component is then removed from the velocity vector and this 
vector is imposed on the boundary. The loss in kinetic energy due to the normal 
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velocity is compensated for by correcting the total energy on the boundary based on 
conservation of total enthalpy. The density gradient, at the wall is assumed to be 
constant. For an // constant wall (constant j), the conservation variables are written 
as 





1 ,A) 

- Q\ 

E'd+l M ^ 



-Q\ 

[ (ij 41 T') ^ 

Q iil.J.k) 

l.fr) 

-Q i 

! (r rJ 41 ,fc) ^ 



1 

2 

r;' 


Vv 

\f>lr + ~i + V: 

x/nl + 'll + vl 

Qlt,j +uk J ! n 


where the magnitude of the normal velocity. l n . is 

Voc , Vv 


^ n u (*j+i,k) 


sfifHif+Ti + r '" J+1M sfvl + + nl ' l '' (,J+1 - k) + n 2 y + n] 

(6.15) 


+ tc. 


(6.14) 




The metric terms ?/ 3 -, ?; y , and r/ 2 are evaluated on the boundary. 


6.4.3 Other 

The following boundary conditions result from computational necessity, not phys- 
ical modeling. 

Changes for Fourth-Order Schemes 

The fourth-order schemes call for special treatment one grid point away from the 
boundaries. Both fourth-order numerical schemes require two points on either side of 
a grid point to compute the derivative at that point. One point from a boundary, the 
scheme must be modified using a skewed difference stencil to prevent over-running 


N AS A/TM— 200 1-210716 


67 


the bounds of the storage arrays and to maintain fourth-order accuracy. The "skewed 
forward" and "skewed backward** fourth-order difference stencils are 


dj _ //+ 3 6 .//+ 2 + 18./,'+ 1 “ 10 /; — *},/?-! 

<Ar /ti-d I2(Ar) 


qi = 3 t /; +1 + io,/;-i8,/;-i + ioA-2-./;-3 

<Ar 12(Ax) 


Pole Boundary 


(6.16a) 


(6.16b) 


The collapsing of a grid surface t o a line or "pole" is oft en done to model cylindrical 
objects in generalized curvilinear coordinates. Figure 6.2 illustrates how a volume in 
computational space is transformed to a sector of a right circular cylinder. Values for 
the variables on the collapsed surface (j m ; n in figure 6.2) are obtained by averaging the 
values over the collapsed index, one point off of the collapsed surface. For example, 
on a constant ?/ surface where ( lines are collapsed to a point, the values of the 
variables at each point on the •// surface are obtained by taking the mass average of 
the variables over all (J at ?/ + 1. This boundary condition is used for the axis in a 
three-dimensional cylindrical calculation. The average flowfjeld values are determined 
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(6.17 


The conservation variables on the boundary are then 



— P avg 



= Pavg 

U avg 


— Pavg 

A’ a vg 

Q 4 [<. J.k) 

— Pavg 

W avg 

Q^ui.k) 

Par 

7 

h + 


. T ,-yPavg (w avg ~f l ’avg "f , 4 avg) 


(6.18) 


Axis of Symmetry 

A boundary condition is needed for axisymmetric analyses on the bounding grid 
line where r = 0, the axis of symmetry. Since the flow is symmetric about this axis 
no flow may cross the grid line. Therefore, the velocity vector must be tangent to the 
boundary. This restriction is the same as the restriction for the slip wall boundary 
and hence, the same method is used to impose the constraints (equation (6.14)). 
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Overlap 


The grid topology commonly referred to as an "(V'-grid is frequently used to 
model round objects and domains, such as cylinders, airfoils and pipes. The do- 
main is formed by wrapping the grid around upon itself so that the minimum and 
maximum boundaries of a common index are coincident (figure 6.3). The coincident 
computational boundaries are in the interior of the physical domain and must be 
treated as if they were also on the interior of the computational domain so that they 
do not create a non-physical barrier on the interior of the domain. In other words, 
the governing equations should be solved on these boundaries as if they are in the 
interior of the computational domain and information must be readily passed across 
the boundaries. This is done by overlapping the grid so that at every point where a 
skewed difference is performed the result can be replaced with a central differenced 
value from a coincident point in the overlapped region. Table 6.1 lists the skewed 
differenced grid points and the coincident central differenced points necessary for the 
overlap boundary condition. 


boundary point 
(skewed difference) 
1 
2 

^max 1 
k max 


interior point 
(central difference) 

h — ^ 

''■mar ’ * 

3 

4 


Table 6.1: Coincident grid points in overlap region 
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6.5 Computer Resource Requirements 


The memory and CPU time required are reported for both the 2D and 3D codes. 
The resource requirements reported are for Silicon Graphics workstations, which were 
used for all the large scale calculations in this document. 

Memory requirements for the two-dimensional/axisymmetric and three-dimensional 
codes are reported in table 6.2. Both single and double precision requitements aie 
given. 


bytes/grid pt. 

code single precision double preci sion 
2D/Axb 244 472 

3D 360 700 


Table 6.2: Memory required 


Table 6.3 contains CPU times for all possible combinations of code (2D, axisym- 
metric, or 3D), equations (Euler, Navier-Stokes, or LES) and numerical scheme. The 
execution times are presented as time per iteration per grid point in seconds. 1 he 
CPU times were obtained using a. Silicon Graphics Power Indigo 2 workstation with 
a MIPS RS000 processor running at 75MHz. The code was compiled using double 
precision data and the compiler options were chosen for maximum performance. 

6.6 Parallel Implementation 

Parallel processing is an efficient and inexpensive method to speed the execution 
of computer programs [63—65]. UNIX workstations with multiple central piocessing 
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rode 

equations 

scheme 

time/iter. /grid pt 

2D 

Euler 

Gottlieb-Turkel 

9.092-10““ 


Euler 

Runge-Kutta 

1.654- 10“ 5 


Navier-Stokes 

Gottlieb-Turkel 

J .218- 10- 5 


Navier-Stokes 

Runge-Kutta 

2. 1 22- J 0 —r ‘ 


LES 

Gottlieb-Turkel 

1 .655- 1 0“ r> 


LES 

Runge-Kutta 

3.196-10“ r ' 

Axi. 

Euler 

Gottlieb-Turkel 

9.4i5-10~ (i 


Euler 

Runge-Kutta 

1.732- 10" 5 


Navier-Stokes 

Gottlieb-Turkel 

1 .498* 1 0 — 5 


Navier-Stokes 

Runge-Kutta 

2.808- 10" 5 


LES 

Gottlieb-Turkel 

1.982- 10“ 5 


LES 

Runge-Kutta 

3.945*10“’’ 

3D 

Euler 

Gottlieb-Turkel 

6.898- 10“ r ’ 


Euler 

Runge-K utta 

1.069-10“'* 


Navier-Stokes 

Gottlieb-Turkel 

1.182-10“ 4 


Navier-Stokes 

Runge-Kutta 

1.555-10“ 4 


LES 

Gottlieb-Turkel 

1 .394- 10 -4 


LES 

Runge-Kutta 

2.074- 10 — 4 


Table 6.3: CPU time required 


units (CPU's) are widely available and cost a fraction of the price of super-computers. 
Codes written to take advantage of the parallel architecture of these machines can 
execute as fast or faster than on a comparable super-computer. Parallel processing 
achieves this performance increase by dividing the work normally done by a single 
central processing unit amongst several CPU’s. 

There are two primary methods used for the parallel processing of CFD solvers. 
The first is distributed memory processing. In this method each CPU uses a separate 
bank of memory for its calculations. Communication of data between CPU's is done 
explicitly through user implemented interfaces. Often the CPU’s do not reside in the 
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same computer, but operate as separate workstations connected through a network. A 
typical distributed memory (’PD solver divides the computat ional domain into seveial 
sub-domains which are each solved on separate processors. The intei faces between 
sub-domains are treated as boundary conditions and are updated by passing updated 
solution values to adjacent sub-domains. This method recpiires limited communica- 
tion between processors, and is therefore suitable for networked workstations, because 
data is only passed at the end of each iteration or set of iterations. It. works very well 
for low order steady-state flow calculations. High order unsteady calculat ions are ill 
suited for t his method because it would be necessary to pass large amounts of data 
between interfaces at every iteration to maintain both temporal and spatial accuracy. 

The second method, shared memory processing, was used in this study. In this 
method, a set of CPU's share access to a common bank of memory. Communication 
between processors is implicit and no special interfaces are required. Shared memory 
systems consist of multiple processors housed in the same computer with a high 
speed high volume bus connecting the processors to the memory bank. Because 
the processors are closely coupled through the shared memory, the division of labor 
between processors can also be closely coupled. This division of labor is done at the 
loop level. The work performed in a. FORTRAN ”D0" loop is divided between CPU's 
with each processor working on a fraction of the total loop. For example, a DO loop 
which is indexed from 1 to 100 could be split between 2 processors with one processor 
operating from 1 to 50 and the second operating from 51 to 100. Implementing 
this method is done through the insertion of compiler directives into the code. The 
directives indicate which loop is to be parallelized and how the variables are shared 
in memory. Care must be taken to identify data dependencies so that variables whose 
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value depends oil the current loop index are not overwritten by processors operating 
on another portion of the loop. Explicit CFD methods are well suited for this type 
of parallelization because the data dependencies are limited and easily overcome. 

The performance of a parallel code is measured by its speedup and efficiency. 
Speedup is simply the ratio of the time required for single processor computation to 
the time required for a multi-processor computation. 

S(«) = E (6.19) 

1 n 


where T n is t he time required for an n processor computation. The efficiency of the 
job is the ratio of the speedup to the number of processors. 


S(w) Ti 

L(n) = = — 

n n 1 ,, 


( 6 . 20 ) 


Ideally a the speedup of a parallel code would vary linearly with the number of 
processors. Two primary factors limit the speedup; processor idle time and commu- 
nication time between processors. Processor idle time is the time the code spends 
outside the parallelized loops. During this time only one processor is active and the 
potential work of the other processors is wasted. Communication time is the de- 
lay during which data is passed between processors and t he shared memory. It is a 
property of the computer and varies greatly between machine types. 

The performance of the 2D/axisymmetric code was measured on a Silicon Graphics 
36 processor Power Challenge computer. The code was run on a representative jet 
calculation using a grid containing 38,829 grid points. A calculation consisting of 
1,000 iterations using the Runge-Kutta scheme was repeated 10 times using an even 
number of processors from 2 to 20. Results were compared to a baseline calculation 
performed on one processor. This procedure was repeated four times. Speedup and 
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efficiency are presented in figures 6.4 and 6.5. The variation in results between each set 
of data is due to the architecture of the computer that was used. The data (memory) 
bus for this machine was insufficiently sized to transfer data for all 36 processors. The 
communication delay time varies with the number of jobs being run on the system 
and the number of processors for a given parallel job. Each set of parallel runs were 
made with the computer under different loads and hence the time required for the 
runs differed. Regardless, the data do show that a linear speedup is seen up to about 
16 processors. Because the computer's performance varied with its load, the parallel 
efficiency cannot be definitively determined. In fact figure 6.5 shows some efficiencies 
exceeding one. This is most likely due to the load on the machine being reduced after 
the one processor job was run. Consequently, the baseline one processor run is not 
appropriate for computing the efficiency of the current multi-processor run. 

6.7 Validation 

To insure accurate implementation of the numerical schemes, the codes were val- 
idated using three simple test cases in both two and three dimensions. The cock s 
results were compared to exact solutions to insure that they produced acceptable 
solutions. 

The standard Runge-Kutta scheme (section 3.2.1) was unstable with or without 
solution filtering and did not provide a converged solut ion for any of the two- or three- 
dimensional calculations. The low dispersion Runge-Kutta scheme (section 3.2.2) 
did provide stable solutions for all calculations attempted. All Runge-Kutta results 
reported hereafter were obtained with the low dispersion scheme (equation (3.9)). 
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6.7.1 Laminar Flat Plate 


Laminar flow over a flat plate was used to check the accuracy of a viscous solution. 
Skin friction along the plate and the self-similar velocity profile are compared to the 
exact solution of Blasius. as given by Schlicting [66]. The skin friction is computed 
using 


c,= 


/' 




( 6 . 21 ) 


and the plate Reynolds number is 


Rtr = 


p ccf oo-r 


( 6 . 22 ) 


The velocity profiles normal to the plate surface can be made self-similar by intro- 
ducing a dimensionless coordinate. 


V = V 



(6.23) 


V vx 

The Reynolds number based on plate length is 10,000 and the freestream Mach 
number is 0.2. 

Two-dimensional solutions were obtained with both numerical schemes. The sixth- 
order filter was used for both calculations. The grid that was used is shown in figure 
6.7. The grid dimensions were 157 points in the streamwise direction and 117 points 
normal to the plate surface. Subsonic inflow 7 and subsonic outflow 7 boundaries were 
specified on the left and right boundaries respectively. A no-slip wall was specified on 
the lower boundary and an extrapolation condition (the same as supersonic outflow) 
was specified on the upper boundary. 

Skin friction results for both the Runge-Kutta and Gottlieb-Turkel scheme are 
shown in figure 6.8. The Runge-Kutta scheme slightly over predicts the skin friction 
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relative to the Blasius solution. The Gottlieb-Turkel scheme underpredicts the skin 
friction and is in greater error than the Runge-Kutta scheme. 

The self-similar velocity profiles are shown in figure 6.9. Several streamwise sta- 
tions. Re r - 2000, 4000, 6000. and 8000 are plotted to verify the code reproduces 
the self-similarity. The Runge-Kutta results are very close to the Blasius solution. 
The Gottlieb-Turkel profiles show a thicker boundary layer thickness with a smaller 
velocity gradient at the wall. This is consistent with the skin friction data. 

To validate the three-dimensional code, a 3D grid was created that consists of 
fifteen evenly spaced 2D grid planes. The 3D code was run using both numerical 
schemes and the results were compared to those of their 2D counterparts. Skin 
friction and velocity profiles at Rt x = 4000 are compared for both schemes (figures 
6.10 T 6.11). The agreement between between the two- and three-dimensional codes 
is excellent for both schemes. 


6.7.2 Supersonic Wedge 


Supersonic flow over a two-dimensional wedge was used to test the ability of the 
code to predict inviscid flows and shock waves. A 15 degree wedge in a Mach 2 
freestream flow was simulated (figure 6.12). Pressure coefficient on the wedge surface 
was compared to the exact solution. 


C P = 


P - Poo 

— f) ( 2 

of'oo'- oo 


(6.24) 


The grid dimensions used were 121 points in the axial direction and bl points in the 
vertical direction. The grid is shown in figure 6.13. Uniform supersonic flow was set 


at the inflow boundary. Supersonic outflow conditions (extrapolation) weie set on 
the right and upper boundaries. A slip wall was specified on the lower boundary. 
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Two-dimensional results are shown in figure 6.14. Overall the results are in good 
agreement with the exact solution. But, the shock location is smeared over several 
grid points and some oscillation in the solution exist both upstream and downstream 
of the shock wave. A comparison of numerical schemes using the sixth-order filter 
shows very little difference. The presence of the shock wave reduces both schemes to 
first-order accuracy (the location of the shock is directly related to the grid spacing) 
nullifying any accuracy advantage of the Runge-Kutta scheme. The affect of fourth, 
sixth, and eighth order filters on the shock location and oscillations was examined 
with the Runge-Kutta scheme (figure 6.14(b)). The lower order filters reduced the 
oscillation without significant affect on shock location or strength. 

Results of the 2D and 3D codes are shown in figure 6.15. For both numerical 
schemes the 2D and 3D pressure distributions compare very well. 

6.7,3 Supersonic Cone 

To test the axisvmmetric terms in the two-dimensional /axisvmmetric code and 
provide a true three-dimensional flowfield for the three-dimensional code, flow over 
a supersonic cone was validated. This test case uses the same flow conditions and 
turning angle as the wedge case. But. the body is treated as a body of revolution 
rather than a planar object. For the axisvmmetric simulation, the grid for the wedge 
case was used (figure 6.13) and the axisvmmetric source terms in the flow solver were 
activated. 

The three-dimensional grid represents a 15 degree section of the cone and its 
flowfield. It was generated by rotating the 2D/axisymmetricgrid about the centerline. 
A 2D/axisymmetric grid plane is located every one degree in the physical domain. 
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The exact solution for the pressure coefficient on the cone surface was obtained 
from the charts in NACA Report 1135, “Kquations, Tables. and Charts lor Com- 
pressible Flow [67]. 

Results for both schemes in 2D are very good (figure 0.17). Both schemes smear 
the shock location over several grid points. Cnlike the wedge calculation ihe shock 
oscillations are weaker and do not over-shoot the pressure behind the shock. The 31) 
results from the Gottlieb-Turkel scheme (figure 6.18(a)) differ from the 2D results 
due to a more pronounced pressure oscillation about the shock. 1 he 3D results from 
the Runge-Kutta scheme (figure 6.18(b)) compare well to its 2D counterpart. 
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(b) internal object in grid (c) grid point specification (filled = ac- 

tive, open = hole, both = boundary) 


Figure 6.1: Modeling an internal object using hole points 
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Figure 6.2: Pol 
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Figure 6.4: Parallel processing speedup 
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Figure 6.5: Parallel processing efficiency 
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Figure 6.6: 2D laminar flat plate 



Figure 6.7: 2D grid for laminar flat plate calculations 
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(a) Gottlieb- Ttirkel scheme 
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(b) Runge-Kutta scheme 

Figure 6.8: Skin friction coefficient for laminar flat plate 
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(a) skin friction coefficient 
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(b) velocity profiles at Rc x — 4000 


Figure 6.10: Comparison of 2D and 3D Gottlieb-Turkel schemes for laminar flat plate 
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(b) velocity profiles at Re x = 4000 

Comparison of 2D and 3D Runge-Kutta schemes for laminar flat plate 
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Figure 6.13: 2D grid for wedge and cone calculations 
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CHAPTER 7 


COMPUTATION OF A NOZZLE FLOWFIELD 

The flow solver described in the previous chapter is now applied to a round super- 
sonic jet. The computational results are compared to high quality detailed data. The 
effects of grid resolution and several numerical modeling parameters on the solution 
are studied. Finally, a thorough examination of the flowfield is presented. 

7.1 Description of the Nozzle 

The nozzle geometry and data of Panda and Seasholtz [68— <1] were selected foi 
study. This data set was obtained using the non-intrusive Rayleigh scattering tech- 
nique. It is based on the measurement of laser light scattered by the air molecules. 
This technique eliminates errors due to probe interference in hot wire and pitot probe 
measurements and biasing errors due to seed particles in Laser Doppler Velocimetiv 
(LDV) and Particle Image Velocimetry (PIV) measurements. The data consists of 
time averaged centerline and radial velocity profiles. 

The nozzle is a one inch diameter convergent-divergent, round nozzle with a de- 
signed exit Mach number of 1.4. The nozzle exhausts into quiescent air. The Reynolds 
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number based oil nozzle diameter is 1.2 million. The nozzle was operated at its min* 
imum shock condition which was slightly less than its design Mach number. The 
operating conditions are given in table 7.1. 


quantity 

symbol 

value 

uni 

ratio of specific heats 

‘■y 

J .4 


nozzle plenum pressure 

IS 

6524.9 

lbs 

ft 2 

nozzle plenum temperat ure 

To, 

540.0 

R 

nozzle exit Mach number 

-V 

1.395 


nozzle exit diameter 

IT, 

0.0833 

ft 

jet velocity 

C, 

1348.4 

s 

ambient pressure 

P 00 

2064.8 

tbs 

ft 2 

an 1 1 ) Sent tern per at u re 


534.6 

R 

Reynolds number 

Re, 

1.2 • 10 6 



Table 7.1: Nozzle operating conditions 


7.2 Flowfield Statistics 

Analysis of RANS solutions is both simplified and limited by the Reynolds aver- 
aging process. The RANS codes readily provide a time averaged solution for analysis. 
However, these solutions lack any unsteady flowfield information. The unsteady solu- 
tions obtained here contain much more information. But the amount of information 
is overwhelming and must be processed and simplified before it can be successfully 
analyzed. 

7.2.1 Time Averaging 

The instantaneous solutions produced by the analysis code sometimes bear little 
resemblance to the time average of the flowfield. To compare the present results 
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with experimental data, empirical correlations and exact solutions a time average is 
computed as 

/(.r) = — ^ f(x.l)dt (<•!) 

and the instantaneous value of the function can be written as a sum of the time 
averaged and fluctuating quantities 

/ = / + /" ( 7 . 2 ) 

In this notation the traditional over-bar ( ) and prime ( ' ) have been replaced with a 
double liar ( = ) and double prime (") t.o avoid confusion with the filtered quantities 
in the LES equations. 

For the jet calculations a time average of density, pressure, arid the three velocity 
components were kept. In addition, the time average of the squares of density and 
velocity components are kept to compute turbulent statistics. 

7.2.2 Turbulent Statistics 

Turbulent statistics can readily be computed from the time averaged flow variables 
being stored during the calculation. 

r f =?-f (7-3) 

By substituting the velocity components into equation (7.3) we obtain the normal 
components of the Reynolds stress u"\ v"\ and w" 2 . The root mean square value is 

simply 
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The root mean square velocities are typically expressed as turbulent intensities. 


Urr. 


a = 


r. u 


v = 


^rrns V 't 




w 


U'rms 

U, 


~w 


t .0 


rx>) 


i . i 


Turbulent kinetic energy is defined as 


1 1 / n2 n 2 \ 

k = 2 1 U + f ’ + «’ J 


(7.8) 


For the LES calculations the turbulent kinetic energy is comprised of resolved and 
sub-grid scale components. 


1 1 / | — " 2 \ , 1 

k — 2 \U + v + w J -|- k sgs 


7.9) 


where 


_ T kk 
' sgs ~ 


7.10) 


We nondimensionalize k using the jet velocity. 


k - = 


i U 2 

2 C J 


(7.11) 


7.2.3 Two Point Correlations 


Correlations of velocity signals at two points within the flowfield can lend insight 
into the structure of the turbulence [72] k [2]. Chu [73] experimentally applied two 
point correlations to jet flows. He obtained a turbulent length scale (eddy size) and 
convection velocity from the correlated data. Scott [74] presented preliminary two 
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point correlation data, obtained with CFD. for Chu's jet. While the comparison of 
the CFD and experimental correlations was inconclusive, due to insufficient simula- 
tion time, this work indicated CFD obtained correlations were were possible. Chu’s 
methodology was followed in this study. 

The location where the turbulent information is desired is specified by the vector 
0i , whose origin is at the center of the nozzle exit (figure 7.1 ). A pair of points evenly 
spaced on either side of this location is specified by a separation Instantaneous 
velocity data at these three points are saved for the correlations. 

The general form of the two point space-time correlation coefficient is 

U'e(o, - V>/2J) Ug [Oi + Vi/2, t + t) 

VC (C,-,t) = = D-i-l 

where Ue is the component of velocity that makes an angle 9 with the jet axis and r is 
a delay time. The correlation coefficient is a measure of similarity of the two signals. 
A value of K near ±1 indicates that the two measurements are highly correlated. 

The location about which the two point correlations were taken was located six 
jet diameters downstream of the nozzle exit on the jet lip line. Eight separations in 
the axial direction of up to one jet diameter were used. Table 7.2 summarizes the 
location and separation vectors for the two point correlations. 

Turbulent Length Scale 

Two point space correlations, where the delay time is zero (r = 0), can be used to 
determine a turbulent length scale. Several sets of correlation coefficients at different 
separations are computed. The separation distance over which the velocities are 
highly correlated is used as an estimate of the length scale L 

(=1 KWi.OWi (7.13) 
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vector coordinates 


local ion 

o, 

[GDj. Dj/2. 0] 

1st separation 

v-; 

[Di/wv 

2nd separation 

*1 

[Dj/7. U. 0] 7 

3rd separation 

'•V 

[^76.0. 0] T 

4th separation 

*7 

[Dj/o.O.Of 

5th separation 

Vi 

[Dj/4,0,0] r 

6th separation 

,M' 

[Djl-l 0,0] r 

7th separation 

Vi 

[I) j /2.Q.0] t 

8th separation 

V? 

[Dj. <U>] r 


Table 7.2: Two point correlation location information 


Convection Velocity 

Two point space-time correlations can be used to determine the convection velocity 
of a turbulent eddy. For a given separation, tpi, the correlation coefficient is computed 
over a range of delay times. The value of r where the correlation coefficient is a 
maximum indicates the time necessary for a disturbance at the upstream point (6 X — 
t/t/2) to reach the downstream point (6 Z + 2). This time and the separation 

distance are then used to compute the convection velocity. 

= ( 7 . 14 ) 

T Um.i* 

7.3 Axisymmetric Solutions 

Assuming t hat the flow is symmetric about, the jet axis reduces the problem to two 
dimensions. Using this assumption would significantly decrease the computational 
expense of a solution and allow for a greater number of parametric investigations to 
be done. However, turbulence is inherently three-dimensional and it is not expected 
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that a solution obtained with the axisvmmetrie version of the ('FT) solver will yield 
physically realistic results. Here, the axisymmetric analyses are used to investigate 
numerical issues, grid resolution, numerical scheme, and boundary conditions, and 
the knowledge gained is applied to the three-dimensional simulations. 

7.3.1 Grid Generation 

The computational grid specifies the nozzle geometry and computational domain 
for the calculations. Unlike previous LES analyses of jets, the internal nozzle contour 
and nozzle lip were modeled. The internal nozzle boundary layer and vortex shedding 
from the nozzle lip may affect the growth and stability of the shear layer. It was felt 
that it was important to include these effects in the simulation. Also, by including 
the internal nozzle the inflow boundary is moved away from the region of inteiest. 
the jet shear layer, reducing the influence of the artificial boundary condition. 

The computational grids used for the analyses were generated using the commer- 
cial software package Gridgen [75]. The grid points were clustered near the nozzle 
walls using hyperbolic tangent stretching to resolve the boundary layers [76]. The 
grid was also clustered near the nozzle exit plane to capture the unsteady vortex 
shedding from the nozzle lip and the initial formation of the shear layer. The cell 
aspect ratio at the upper and lower corners of the nozzle hp was set to one. Grid 
point clustering in the expected region of the shear layer downstream of the nozzle 
lip was also implemented. The computational domain extends 20 nozzle diameters 
downstream of the nozzle exit in the axial direction and 10 nozzle diameters from the 
jet centerline in the radial direction. A representative grid with 301 point in the axial 
direction and 129 points in the radial direction is shown in figure 7.2. 
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7.3.2 Specification of Boundary Conditions 

All the axisy linnet ric jet computations used the same boundary conditions (figure 
Subsonic inflow using the jet plenum conditions was specified at the inflow of 
the internal nozzle section. The no-slip wall condition was used for the internal and 
external nozzle surfaces. An axis of symmetry/ slip wall condition was used on the jet 
centerline. 

Simulating the quiescent conditions of the surroundings can cause numerical dif- 
ficulties. 1 he boundary conditions are formulated assuming a known flow direction 
(inflow or outflow). In the still air small disturbances may cause portions of a bound- 
ary to have an inflow and other portions to have an outflow. This leads to over or 
under specification of the boundary and numerical errors result. To overcome this 
problem a small freestream flow is imparted to ensure a coherent flow direction and 
properly specified boundaries. In this case a Mach 0.05 freestream was used. Subsonic 
inflow was specified on the upstream external boundary using the freestream total 
conditions. Conditions on the upper boundary were extrapolated from the freestream. 

The subsonic outflow boundary condition imposes a constant pressure over the 
entire boundary. In reality the pressure on the outflow boundary varies in both time 
and space. To accommodate this pressure variation while maintaining the correct 
freestream flow a combination of subsonic outflow and extrapolation conditions were 
used. The subsonic outflow condition was specified on the upper portion of the down- 
stream boundary to maintain the correct pressure level in the freestream. Conditions 
on the rest of the downstream boundary were extrapolated from the interior to ac- 
commodate both temporal and spatial pressure variations near the jet centerline. 
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7.3. 3 Effect of Grid Resolution 


The effect of grid resolution on the solution was investigated. Three grids were 
used, where the grid spacing in each direction was halved from the previous grid. All 
solutions were obtained with the Runge-Kutta scheme using an eighth-order filter. 
The sub-grid scale model was not used. Because the grid spacing is not uniform, 
it is difficult to quantify. Here, the spacing was represented by a nondimensional 
computational grid size. 

Aa'* = • At/ (7-15) 

As the grid resolution was increased, the resolution of the unsteady flowfield im- 


grid 

dimensions 

A.r* 

Umar 

I'max 

k* 

■m a.r 

coarse 

151 x 

65 

CN 

O 

o 

O 

0.20006 

0.15480 

0.047999 

medium 

301 x 

129 

5.1031 • 10-' 3 

0.26266 

0.20305 

0.096392 

fine 

601 x 

257 

2.5516 • 10-' 3 

0.32909 

0.22788 

0.12178 


Table 7.3: Effect of grid resolution on axisymmetric solution 


proved. The change in entropy of the flowfield can be used to visualize the tui Indent 
structures in the flowfield. The viscous mixing of the jet and ambient air increases the 
entropy of the flow. The vortices alter the shape of the mixing region and this altered 
shape can be seen in the gradients of entropy in the mixing layer. Figure 7.4 shows 
contours of entropy for the three grids used. The plots show a dramatic increase in the 
resolution of vortical structures with grid resolution. A more quantitative measure, 
the maximum of the turbulent statistics in the flowfield. are shown in figure 7.5. The 
turbulent intensities and kinetic energy increase as the grid spacing decreases (as the 
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grid 

dimensions 

a7-/k 

coarse 

151 : 

c 65 

2081.5 

medium 

301 j 

< 129 

1040.8 

fine 

601 : 

c 257 

520.40 


Table 7.4: Comparison of grid spacing to the Kolmogorov scale 


grid resolut ion increases). This result is expected as the intensities and kinetic energy 
should increase as the grid is refined and the scheme increasingly captures more of 
the smaller energy containing eddies. With even further resolution one could expect 
the energy to decrease as the inertial subrange and Kolmogorov scales are reached, 
where the energy is dissipated. 

To compare the grid spacing to the Kolmogorov scale, an average cell size was 
computed as follows 

\J Lgrid ' H grid ( t .16) 

where L grt( ] and H grtt { are the overall length and height of the grid. The Kolmogorov 
scale was estimated using equation (1.1) and assuming the integral length scale is 
approximately \D 3 , Table 7.4 shows that even for the finest grid the cell size is over 
500 times larger than the Kolmogorov scale. 

7.3.4 Comparison of Numerical Schemes 

A comparison of the two numerical schemes was performed on the medium sized 
grid (301 x 129). The solutions were run for two characteristic acoustic times (the 
time for an acoustic wave to pass through the solution domain) to establish proper 
initial conditions for obtaining turbulent statistics. The solutions were then run for an 
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additional two characteristic times and the flow field was averaged. Both schemes were 
run using their maximum stable time step. Figure 7.6 shows the entropy contours for 
the two schemes. The Runge-Kutta scheme exhibits more vortical structure and the 
initial vortex roll-up of the shear layer occurs earlier. Table 7.5 quantifies the results in 
terms of turbulent intensities and kinetic energy. The Runge-Kutta scheme predicts 
higher maximum turbulent intensities and kinetic energy than the Gottlieb- Turkel 
scheme, indicating better resolut ion of the flowfield. I lie simulation time presented 
in the table is the ratio of the time required to the time required for the Gottlieb- 
Turkel scheme. Consistent with the one-dimensional results, the Runge-Kutta scheme 
is more computationally efficient, using 16.5 percent less ( PI time. 

scheme Umax t'mru- k m <ij- time 

Gottlieb- Turkel 0.26033 0.18561 0.080193 1.00000 

Runge-Kutta 0.26266 0.20305 0.006302 0.83523 


Table 7.5: Effect of numerical scheme on axisymmetric solution 


7.3.5 Evaluation of Exit Zone Boundary Condition 

The use of an exit zone outflow boundary has been advocated in the area of 
computational aeroacoustics where proper resolution of very low magnitude acoustic 
waves is critical. It is not clear that an exit zone is required when the simulation is 
not intended to capture acoustic waves. The need for and effectiveness of the exit 
zone boundary condition was tested by comparing solutions with and without this 
boundary treatment. The computational grid was modified by adding an exit zone. 
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a region of highly stretched grid cells downstream of the previous outflow boundary 
(figure 7.7). Fifteen additional grid planes were added using a ten percent geometric 
stretching factor. The simulations were run using the Runge-Kutta scheme with a 
sixth-order filter. 

Entropy contours show some differences in the turbulent structures. Blit , it is not 
c lear that the differences are due to any reflecting waves. Turbulence levels given in 
table 7,6 indicate very little difference between the two solutions. 

boundary condition u max v max k; ?KU . 

Exit Zone 0.30181 0.20379 0.096595 

No Exit Zone 0.26266 0.20305 0.096392 

Table 7.6: Effect of exit zone boundary condition on axisymniet ric solution 


7.3.6 Effect of the Sub-Grid Scale Model 

The solutions obtained thus far have not used the sub-grid scale model. They are 
in effect DXS solutions. However, the grids used are not fine enough to resolve all 
the turbulent scales. Scales smaller than the grid size are not computed and their 
contribution is lost. Since these scales tend to dissipate the larger eddies, the solu- 
tions without the sub-grid model should predict larger more energetic eddies than a 
corresponding LES solution. The use of the term DNS to describe these solutions is 
misleading because DNS implies that all the turbulent scales are computed. There- 
fore, these solutions will be referred to as “coarse grid DNS" solutions. 
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An LES solution on the medium grid was run and compared to its “coarse grid 
DNS" counterpart. Entropy contours are compared in figure 7.9. The additional 
dissipation from the sub-gricl model serves to damp the large scale structures. This 
behavior is evident by observing the initial vortex roll-up in the mixing layer. The first 
large eddies occur further downst ream in t he LES solution due to the eddy viscosity 
from the sub-grid model. Overall, the LES solution produces less turbulent mixing. 
7 percent less turbulent intensity and 14 percent less turbulent kinetic energy (table 
7.7). 


small scale modeling u max ^ ma.r 

without sub-grid model 0.26266 0.20305 0.096392 

with sub-grid model 0.24247 0.18797 0.082770 


Table 7.7: Effect of sub-grid scale model on axisymmetric solution 


7.3.7 Evaluation of the Axisymmetric LES Solution 

To obtain a steady time averaged solution, the LES solution was run for 100,000 
iterations, which corresponds to 0.0052 seconds of simulation time. The instantaneous 
solution was averaged every 100 iterations. 

Analysis of the Flowfield 

Contour plots of sample instantaneous and time averaged flowfield quantities are 
shown in figures 7.10 through 7.12. The turbulent structures present in the mixing 
layer are clearly seen in the instantaneous contour plots. The density' contours (figiuc 
7.10) reveal the shock structure in the potential core. The radial velocity contours 
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(figure 7.1‘2) provide good visualization of the vortical structures. The time averaged 
plots show smoothly varying gradients in the mixing layer and are analogous to a 
flowfield obtained from a RANS solution or analytic method. 

Plots of p. u. v. and k. are show in figures 7.13 through 7.16. All the turbulent 
quantities exhibit a similar structure. The detail plots of the nozzle exit show a 
small concentrated unsteady region at the nozzle lip caused by vortex shedding. This 
region is followed by a region dominated by very small scale turbulence. Radial 
velocity contours in this region (figure 7.17(a)) show that no large scale structures 
exist. But. a plot of the sub-grid scale turbulent kinetic energy, t**, shows that the 
energy from the small scales peak in this region (figure 7.17(b)). Downstream of the 
small scale region, the mixing layer becomes unstable and large vortical structures 
begin to form and grow larger with increasing distance from the nozzle exit. 

Comparison to Experimental Data 

fhe time averaged solution is compared to the experimental data in figures 7.18 
and 7.19. The experiment al data exhibits the behavior of a. typical jet [77,78]. A series 
of weak shock waves can be seen near the nozzle exit. The end of the potential core, 
as indicated by the start of the velocity decay on the centerline, is at approximately 
7.5 jet diameters. The predicted velocity on the jet centerline captures the first 
few shock waves near the exit. The remaining shocks are not captured due to grid 
stretching away from the nozzle exit. The predicted centerline velocity does not decay 
within the computational domain indicating that the axisymmetric assumption has 
constrained the mixing layer to lie above the jet axis. Radial profiles of axial velocity 
at :r / D } — 2. 4. 6. 8. 10. and 12 are shown in figure 7.19. The experimental data 
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was taken both above and below the jet axis. All these data are plotted to show the 
asymmetry in the data. For comparison the LES solution is reflected about the axis. 
Agreement between the CFD and experiment is very good up to G jet diameters. 
Beyond this point the prediction departs from the experiment due to the lack of 
decay of the potential core. The peak in the experimental profiles begins to decrease 
and the jet spreads at a greater rate than the CFD indicates. The good agreement 
upstream of the end of the potent ial core indicates that the axisymmetric assumption 
is valid in this region. In the region where the potential core breaks down, highly 
three-dimensional turbulent structures may be the cause. 

7.4 Three-Dimensional Solutions 

A large-edcly' simulation of the full three-dimensional jet flowfield is a large un- 
dertaking. The number of grid points and time necessary for a solution severely limit 
the size and number of calculations. The results from the axisymmetiu solutions, 
code validation, and one-dimensional error analysis were used to guide the simulation 
process. The Runge-Kutta scheme has been proven superior to the Gottlieb-1 urkel 
scheme and was used for all 3D computations. 

Three computations were performed. The first used the sixth-ordei filtei and 
provided time averaged flowfield data. The second calculation used the eighth-order 
filter and was used to ascertain the effect of the solution filter and to obtain turbulent 
statistics and a series of instantaneous velocity distributions. The third used the 
sixth-order filter and a modified grid to obtain two point correlation information. 
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7.4.1 Grid Generation 


The axisvmmet lie grid study showed an increase in resolution of the turbulent 
structures with increasing grid resolution. But. a trade between accuracy and com- 
putational cost must be made to keep the simulation time reasonable. The resolution 
used in the fine grid would create a prohibitively large 3D grid and the coarse grid 
showed no visible turbulent structures. The medium grid (301 x 129) showed reso- 
lution of large scale turbulence and provides a basis for a reasonably sized 3D grid. 
Also, a reasonable sub-grid model should replace the effect of the unresolved terms. 
To form the three-dimensional grid, a grid plane corresponding to the axisymmetric 
grid was positioned at every ten degrees around the jet axis resulting in a cylindrical 
domain. Four additional grid planes were added overlapping the first four planes 
to facilitate use of the overlap boundary condition. The final grid has dimensions 
301 x 129 x 40 and contains 1.553.160 points. Figure 7.20 show r s streamwise and 
cross-stream grid planes. 

The grid was modified slightly for the two point correlation work. The grid points 
were redistributed locally so that a grid point was located at each correlation location 
(table 7.2). The overall grid structure and size remained the same. 

7.4.2 Specification of Boundary Conditions 

Boundary conditions are specified in the same manner as done in the axisymmetric 
case (section 7.3.2). No-slip walls are specified on the nozzle surfaces. Subsonic inflow 
conditions are used for the nozzle plenum and freestream inflow boundaries. The 
combination of extrapolation and subsonic outflow' are used on the outflow plane to 
allow for pressure variation and extrapolation is used on the upper boundary. 
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Two additional boundary conditions are necessary for the 3D grid. First the 
overlap condition is used to create a continuous domain in the azimuthal direction. 
Second, the jet centerline consists of a grid plane that is collapsed to a line. The pole 
boundary is used here. 

7.4.3 Presentation of the Time Averaged Flowfield 

The LES simulation was run until the time averaged centerline velocity profile 
remained unchanged for over 5.000 iterations. This required 50.500 iterations using 
an average time step of approximately 48 • 10 _y seconds. Using 10 processors on the 
Silicon Graphics Power Challenge machine, the simulation required about two months 
of calendar time for completion. The solution was averaged every 100 iterat ions. 

Contours comparing instantaneous and time averaged velocit ies in both the stieam- 
wise and cross-stream planes are presented in figures 7.21 - 7.26. 1 he difference be- 
tween the instantaneous and time average is striking. The instantaneous velocities 
show large turbulent structures that, average to zero over time. Large azimuthal ve- 
locities and the appearance of turbulent structures that cross the jet axis near the 
end of the potential core indicate a highly three-dimensional flowfield. 

Radial velocity and sub-grid turbulent kinetic energy near the jet lip are shown 
in figure 7.27 and differ from the axisymmetric solution (figure 7.17). In the mixing 
layer the vortical structures in the 3D solution appear weaker because of the 3D relief 
effect not accounted for in the axisymmetric solution. Also, the 3D solution does not 
show any resolved vortices near the jet lip. The sub-grid turbulent kinetic is much 
higher, in the 3D solution, indicating that the sub-grid model has dissipated this 
motion. The increased effect of the sub-grid model may be caused by the large grid 
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sparing in the azimuthal direction. This large grid spacing results in a larger filter 
width, which increases the amount of eddy viscosity (equation (5.19)). 

Contours of dilatation. 0S are shown in figure 7.28. By rewriting the continuity 
equation (equation 2.1) we can see that dilatation is related to the convection of 
density waves, which can be related to the convection of sound waves. 


dit, _ 1 (dp_ dp \ 

dx; P V dt Ul d-Xi J 


(7.17) 


The dilatation contours show that the sound for this jet would emanate from the 
mixing layer near and the end of the potential core. 

The time averaged velocities closely resemble both experimental data and RANS 
calculations. But. clearly the LES solution is capable of providing much more insight 
into the flow physics of the jet through the unsteady information. 

The time averaged velocity profiles on the jet centerline are presented in figure 
7.29. The LES solution predicts that the potential core is shorter than found experi- 
mentally indicating that the turbulent eddies are too energetic. The sub-grid model 
may not be providing adequate dissipation of the large scale eddies. Adjusting the 
constants to the model may correct this problem. The rate of velocity decay beyond 
the potential core is very close to that of the experimental data. 

Radial profiles of axial velocity at x/Dj = 2. 4. 6. 8, 10, and 12 are shown in 
figure 7.4.5. At x/Dj — 2 and 4 the 3D prediction shows less spreading than either 
the experimental data or the axisymmetric prediction. This difference may be caused 
by the truncation error due to the large grid spacing in the azimuthal direction. 
The spacing in the azimuthal direction is much larger than the spacing in the axial 
and radial directions and is the reason for the difference between the axisymmetric 
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and 3D solutions. Since the grid resolution is not sufficient to capture the turbulent 
structures in this region, the sub-grid model should provide the equivalent effect 
through increased eddy viscosity. The fact that the LES solution predicts less jet 
spreading than the experiment indicates that the sub-grid model is not adequate. 

While the agreement with the experimental data is far from perfect, this calcu- 
lation is superior to the other high Reynolds number LES calculations in the litera- 
ture [31-33]. 

7.4.4 Presentation of the Instantaneous Flowfield and Tur- 
bulent Statistics 

A second nearly identical simulation was run to obtain turbulent statistics and in- 
stantaneous flowfield data not saved during the first simulation. The order of the filter 
was changed from sixth to eighth to ascertain its effect. The simulation was again run 
until the centerline velocity profile did not change over a period of 5.000 iterations. 
The time step was fixed at 50- 10“ 0 seconds and 70.200 iterations were required. The 
solution was averaged every 80 iterations. The increase in the filter s order reduced 
the amount of numerical dissipation in the solution and increased the resolution of 
the high w r ave number disturbances. This leads to increased turbulent mixing as 
shown in the centerline velocity profile (figure 7.31). The increased turbulent struc- 
tures are most likely the reason why more simulation time was required to obtain a 
suitable average. In theory the sub-grid scale model should adjust to the change in 
flowfield resolution and filter and grid independence should be possible. The increase 
in resolution of the turbulent eddies, with a corresponding increase in |.Sj should be 
compensated for by an increase in eddy viscosity (equation 5.19). It is evident that 
the current model and the choice of coefficients are not adequate for this problem. 
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Instantaneous velocity contours at three times are shown in figures 7.32 - 7.37. 
The contours indicate that the jet flowfielcl is highly three dimensional with very 
large scale structures. Significant variation with time is found not only in the mixing 
layer, but also in the potential core. Here the core flow "flaps" and rotates about 
the jet axis. Also, the shock structure varies considerably due to the variation in the 
boundary between the core flow and mixing layer caused by the vortical motion in 
the mixing layer. 

Turbulent intensities and turbulent kinetic energy are presented in figures 7.38 
and 7.39. The axial turbulent intensity peaks at u = 0.4123 in the mixing layer 
near the nozzle lip. Unsteady motion is also found in the potential core, due to 
the unsteady shock structure, and in the freest-ream, due to acoustic waxes. The 
radial and azimuthal intensity plots are lower in magnitude and peak at a value of 
ft = 0.2002 at the end of the potent ial core. They have very similar structures except 
for two very low intensity regions. The first, is the unsteady radial component of 
velocity induced by the shocks in the potential core. The second is a wide region of 
unsteadiness in radial velocity at the downstream boundary. 

Turbulent kinetic energy is plotted in figure 7.39. The peak values, k* = 0.1760. 
occur in the mixing layer just downstream of the nozzle lip. A second lower peak, 
k* = 0.1210. occurs near the end of the potential core. 

The anisotropy of the turbulence is easily seen in figure 7.40. The ratio of radial 
to axial turbulent intensity, u/v, is shown in three plots with different contour levels 
to separate the regions of isotropy and anisotropy. The vast majority of turbulent 
structures are anisotropic with a peak ratio of 2.775. The majority of the fiowfield 
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has a ratio in the range J.2 < u/v < -.0 1 his result indicates that the isotropic 
assumptions of most RANS turbulence models are not applicable for this flow. 

7.4.5 Presentation of the Two Point Correlation Data 

A simulation on the modified 3D grid was run using the sixth-order filter and axial 
and radial velocity data were sawed for correlation. These velocity components can 
be combined to form the velocity vector making an angle theta with the jet axis. 

( >g — a cos 0 + r sin 6 ( 7 . 18) 

A time history of the velocity signals at the points for the fourth separation, \/2D ,, 
arc* shown in figure 7.4.5. Three angles are examined 0. 45, and 90 degrees. There 
is very little change in the velocity signal as it convects from the upstream to the 
downstream point. The average, velocity and turbulent intensity for each angle is 
given in table 7.8 

9 Oe/Uj 

0 0.65417 0.26473 

45 0.46029 0.26601 

90 -0.0032162 0.17076 

Table 7.8: Velocity statistics at the two point correlation location 


Two Point Space Correlation 

A plot of the two point correlation coefficients versus separation distance is shown 
in figure 7.42. The correlation of the velocity signal decreases towards zero with 
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increasing separation distance indicating that the average size of the turbulent struc- 
tures is less than the maximum separation distance (one jet diameter). The turbulent 
length scab 4 was estimated using equation (7.13). To evaluate the integral, the data 
was fit with a sixth-order polynomial using the method of least squares [79]. The 
polynomial is of the form 

y = C(s.r h -f C 5.?' 0 + C 4 X { + c* 3 .r ^ + c* 2 .r“ + c\.v + cq ( / .19) 

The coefficients obtained in the curve fitting are found in table 7.9. The polynomial 


coefficient 


angle 

<70 

t'l 


c?. 

C.3 

1*5 

CG 

h norm 

0 

1 . 0 UU 

0.041 24 

-4.230 

5.458 

-1.940 

-1.458 

1.105 

8.564 -10” 1 ' 

15 

1 .004 

-0.01181 

-5.249 

6.480 

-1.736 

- 2.002 

1.305 

9.989 U0“ p 

90 

1 .003 

0.03630 

■11.13 

1 6 .46 

-2.633 

-8.769 

4.687 

4.052 10 ~ 5 


Table 7.9: Curve fit coefficients for two point space correlations 


was then integrated analytically to obtain the length scale. The results are shown in 
table 7.10. The data at 0 degrees indicate that the extent of turbulent structures in 


0 ( Vh 

0 0.50161 

45 0.37032 

90 0.10422 


Table 7.10: Turbulent length scales 


this direction are approximately a half of one jet diameter. The smaller length scales 
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obtained from the 45 and 90 degree data are caused by the relatively large region of 
negative correlation coefficients for this data. This effect was also noted by Chu [74]. 

Two Point Space-Time Correlations 

Plots of the two point space-time correlations for each separation are shown in 
figure 7.4.5. For each separation there is a delay time where the signals are clearly 
correlated (a distinct peak in the correlation curve). This delay time increases with 
increasing separation distance. 

The exact location of the maximum correlation coefficient cannot be determined 
from the discrete data obtained in the analysis. To find this location, the correlation 
data was again curve fit with a polynomial (equation (7.19)). The maximum in 
the curve is found by differentiating the polynomial equation. The location of the 
pea.k in the correlation curves corresponds to the time required for a. signal at the 
upstream point to con vec t to the downstream point. The resulting convection speed. 
U c , (equation (7.14)) is shown in figure 7.44. The convection speed decreases slightly 
with increasing separation distance due to viscous effects. An average of the computed 
convection speeds is 0.62607 Jj which is consistent with the results of Chu [73]. 
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Table 7.11: Curve fit coefficients for two point space-time correlations 
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Figure 7.2: Axisvmmetric computational grid 
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Figure 7.3: 


Boundary conditions for axisymmetric 


calculations 
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Figure 7.5: Turbule 
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Figure 7.7: Axisymmetric computational grid with exit zone 
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(b) without exit zone 


Figure 7.S: Entropy contours for exit zone boundary condition comparison 
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(b) with sub-grid model 


Figure 7.9: Entropy contours showing effect of the sub-grid model 
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(b) time averaged 


Figure 7.10: Density contours for the axisymmetric solution 



(a) instantaneous 



(b) time averaged 


Figure 7.11: Axial velocity contours for the axisymmetric solution 
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(b) time averaged 


Figure 7.12: Radial velocity contours for the axisymmetric solution 



(b) detail of nozzle exit 


Figure 7.13: Root mean square density contours for the axisymmetric solution 
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(a) entire domain 



(b) detail of nozzle exit 


Figure 7.14: Root mean square axial velocity contours for the axisymmetric solution 



(b) detail of nozzle exit 


Figure 7.15: Root mean square radial velocity contours for the axisymmetric solution 
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(a) entire domain 



(b) detail of nozzle exit 


Figure 7.16: Turbulent kinetic energy contours for the axisymmetric solution 
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Figure 7.18: Time averaged centerline velocity profile for the axisymmetric LES 
lution 
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(a) x/Dj = 2 
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(b) x/Dj = 4 


Figure 7.19: Time averaged radial profiles of axial velocity for the axisvmmetric LES 
solution 


continued 
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(a) streamwise, x-y plane 



(b) crossstream. v-z plane at outflow boundary 


Figure 7.20: Three-dimensional grid 
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(a) instantaneous 



(b) time averaged 


Figure 7.21: Axial velocity contours for 3D LES solution 


N AS A/TM— 200 1-210716 


136 



(b) time averaged 


Figure 7.22: Radial velocity contours for 3D LES solution 
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(a) instantaneous 



(b) time averaged 


Figure 7.23: Azimuthal velocity contours for 3D LES solution 
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(a) instantaneous 


(b) time averaged 


Figure 7.24: Total velocity contours at x/Dj = 3 for 3D LES solut ion 



Figure 7.25: Total velocity contours at x/Dj = 6 for 3D LES solution 
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Figure 7.26: Total velocity contours at x/D } 


9 for 3D LES solution 
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(a) radial velocity 



(b) sub-grid turbulent: kinetic energy 


Figure 7.27: Nozzle exit detail for 3D LES solution 
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Figure 7.32: Instantaneous axial velocity contours for the 31) LES solution 
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(c) jff = 32.297 


Figure 7.33: Instantaneous radial velocity contours for the 31) LES solution 
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Figure 7.34: Instantaneous azimuthal velocity contours for the 3D LFS solution 
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(a) = 0.0000 
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(b) 



16.149 



32.297 


Figure 7.35: Instantaneous total velocity contours at x/Dj = 3 for the 3D LES 
solution 
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Figure 7.36: Instantaneous total velocity contours at x/Dj = 6 for the 3D LES 
solution 
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Figure 7.3' 
solution 



(a) ^ = 0.0000 
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(b) radial turbulent intensity 



(c) azimuthal turbulent intensity 


Figure 7.38: Turbulent intensity contours for the 31) LFS solution 
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Figure 7.39: Turbulent kinetic energy contours for the 3D LES solution 
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(a) u/v < 0.8 



(b) 0.8 < S/r < 1.2 



(c) u/t? > 1.2 


Figure 7.40: Ratio of axial to radial turbulent intensity for the 3D LES solution 
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Figure 7.41: continued 



ax/d. 

Figure 7.42: Two point spB.ce correction coefficient 
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(b) 45 degrees 


Figure 7.43: Two point space-time correlation coefficient 


continued 
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Figure 7.43: continued 



Ax/D. 


Figure 7.44: Convection velocity 
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CHAPTER 8 


CONCLUSIONS AND RECOMMENDATIONS 


A numerical method to simulate high Reynolds number jet flows was formulated 
and applied to gain a better understanding of the flow physics. Close attention is 
paid to the sources of error in such calculations and efforts were made to minimize 
them whenever possible. 

Large-eddy simulation is chosen as the most promising approach to model the 
turbulent structures due to its compromise between accuracy and computational ex- 
pense. The filtered Xavier-Stokes equations are developed including a total energy 
form of the energy equation. Sub-grid scale models are adapted from compressible 
forms of Smagorinsky's original model. 

The effect of using disparate temporal and spatial accuracy in a numerical scheme 
was discovered through one-dimensional model problems. The lower order time step- 
ping found in many schemes, such as the Gottlieb- Turkel scheme examined here, 
causes the scheme to revert to the lowest order accuracy. A new uniformly fourth- 
order accurate numerical method was developed based on this woik. The scheme 
consists of a low-dispersion Runge-Kutta time stepping scheme with a central diffei- 
enc.e spatial operator. Solution filtering is used to maintain stability. Results in both 
one- and two-dimensions showed this new scheme clearly superior to the second-order 
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iii time and fourth-order in spare Gottlieb- Turkel scheme. The measure used to judge 
t he schemes was the computational efficiency, the time required to reach a given level 
of error. 

1 he resulting flow solver was configured to run on a shared memory parallel com- 
puter. Poor computer architecture prohibited conclusive results, but the evidence 
indicates that the code scales well up to 16 processors. Results from validation ex- 
ercises show that the code accurately reproduces both viscous (laminar flat plate) 
and inviscid (supersonic wedge and cone flows) flows. The validation exercises also 
confirmed the increased accuracy of the new numerical scheme. 

Numerous axisymmetric simulations were performed to investigate the effect of 
grid resolution, numerical scheme, exit boundary conditions and sub-grid scale model 
on the solution. While the axisymmetric assumption was not accurate for the jet 
flowfield, valuable information was gained for use in the three-dimensional caleula- 
t ions. 

The three-dimensional calculations showed that this LES simulation accurately 
captures the physics of the turbulent jet. The agreement with experimental data 
relatively is good and is much better than results in the current literature. However, 
there is still much room for improvement. The improved agreement over previous 
work can be attributed to the new numerical scheme and the modeling of the nozzle 
lip. A subsequent run using a higher-order solution filter indicated that the modeling 
of the unresolved scales needs improvement. 

Several techniques were used to gain a better understanding of the underlying 
physics. A plot of dilatation was used to provide insight into the sound field by indi- 
cating the location of acoustic sources in the jet mixing layer. Turbulent intensities 
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indicate that the turbulent structures at this level of modeling are not isotropic and 
this information could lend itself to the development of improved sub-grid scale mod- 
els for LES and turbulence models for KAN'S simulations. A two point correlation 
technique was used to quantify the turbulent structures. Two point space correla- 
tions were used to obtain a measure of the integral length scale, which proved to 
be approximately \D 3 . Two point space-time correlations were used to obtain the 
convection velocity for the turbulent structures. This velocity ranged from 0.57 to 
0.71 Vj. 

There are several recommendations for further work. 

The accuracy of the simulations is highly dependent on grid resolution. Accurate 
resolution of the shock structure in the potential core was found with the axisymmetric 
calculation. However, the large grid spacing in the azimuthal direction in the ID 
calculations diminished this. Further resolution in this direction may improve the 
prediction of both the shock structure prediction and the mixing layer. A more 
systematic study of grid resolution in all three directions is desirable, but may be 
computationally prohibit i v e . 

A change in the order of the solution filter drastically changed the centerline ve- 
locity decay. The change in turbulent mixing with an increase in resolution of scales 
indicates that the sub-grid scale model is not accurately mimicking the effects of the 
unresolved scales. Further research into improving the sub-grid scale models is neces- 
sary. A promising approach is the dynamic sub-grid model [26] which automatically 

adjusts the coefficients based on the filter width. 

The time required to run a simulation on this relatively simple geometry severely 
limits the usefulness of LES. This limit may be eased if a more efficient numerical 
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method was found. As shown, the time stepping scheme is typically the factor that 
limits the computational efficiency. A high-order accuracy efficient time stepping 
scheme would allow faster turn-around of solutions and more accurate answers on 
finer grids. 

While there is room for improvement in accuracy, this research has shown that 
large-eddy simulation can be used, as is. to provide new insight and information about 
high Reynolds number jet flows. The characterization of the turbulent structures, size, 
convection speed, and degree of anisotropy, can be used to develop improved tools 
for predicting the fluid mechanics and acoustics of jets. 
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APPENDIX A 


DERIVATION OF THE FILTERED EQUATIONS 

A complete derivation of the Favre filtered Navier-Stokes equations used in Chap- 
ter 5 is presented. 

A.l The Filter 

A filtering function 6', is used to separate large and small scale components. The 
filtering operation applied to a function / is 

/ = G(x - 0/(0<# 

The function can then be decomposed into its resolved /filtered, /, and unresolved, / 
parts 

/ = / + /' ( A - 2 > 

For most applications, the function is not specified. But, several constraints are 
placed on the the function to ensure that the filter commutes with the derivative. 

W. = ( 2l (A. 3) 

dx dx 
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The constraints are 


1) &(-{) = G(£) 

2) / ' G(0<1£ = I 

3) G'(£) -» 0 as |£| — * oc 

4) G(£) is small outside ( — y. y) 

where A is a characteristic width of the filter function. Where it has been necessary 
to know the form of the filter function, researchers have typically used either a box. 
Gaussian, or spectral cutoff filter. 


Favre (density) weighting is used in the filtering process. This allows for convenient 
recovery of terms corresponding to the unfiltered equations. 



(A. 4) 


A. 2 Continuity Equation 


F iltering the continuity equation is a straight forward process. The spatial filter 
is first applied to the continuity equation (2.1). 


dp dpu, 

— 4- -2—2 = 0 

8i ^ dxi 

Since the filter commutes with the derivative equation (A.o) is rewritt 


(A.5) 


en as 


dp dpu • 
dt dx. 


(A. 6) 


Then F a vre weighting is used to recover an equation of the same form as the unfiltered 


equation (2.1). 


dp dpu i _ 

di dxi 


(AS) 
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A. 3 Momentum Equation 


The filtering of the momentum equation and definition of the sub giid scale stress 
tensor is presented below. The filtering operation is applied to the momentum equa- 
tion (2.2). 

dpy, tipUjVj dp _ (A. 8) 

dt c)xj d.v, dx :i 

Using the the property in (A. 3) the equation is rewritten 

(fpy t d puj-u j dp _ d&ij (A. 9) 

dt dxj dxi dxj 

Favre weighting is then applied 

dpu, dpujn, dp _ dajj ( A. 10) 

dt + dxj dxi dx.j 

where the filtered stress tensor is 

<t ij - — § pSijSkk + 2 pS tJ (A.ll) 

The filtered stress tensor and the term upij are not in useable forms because they are 
the filter of a product of two variables. We define a new resolved stress tensor as 

= ~lp^S kk + 2/1 Su- ( A - 12 ) 


where jj = n{T) and 


• _ 2 V c)x t Oxj J 


(A. 13) 


We can then write (ATI) as 


= Vij + W13 


(A. 14) 
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Similarly, we writ* 


piliUj = pV,Uj + /) ( 11,11, - u t u ■) 


We define the sub-grid scale stress tensor as 


• j = P ( 


u,ii , — u, u. 


and rewrite (A. 15) as 


pupij = pu,iij 4- T, 


(A. 17) 


Finally, we substitute (A. 14) and (A. 17) into (A. 10) and obtain the final form of the 
filtered momentum equation 


dfm, t dpUjUj c)p d(Tij dra i d ^ 

~dF + dxj + ^ = ~d7~ ~ + {(T " ~ a ‘-' : 


(A. 18) 


A. 4 Energy Equation 

The filtering the energy equation is the most involved process. Because the total 
energy form of the energy equation is used, several additional manipulations of the 
resulting terms are required in order to recover terms for which there are sub-grid 
scale models. Most work in compressible large-eddy simulations have used either 
static or total enthalpy forms of the energy equation. First the filtering operation is 
applied to the energy equation (2.6) 


dpe t dpipe, du t p _ du 3 i 7 tJ dq , 

di dx, dxj dx t dx t 

Commuting the filter operation with the derivative we obtain 

dWt , dWdp duiP _ _ dfh_ 

dt dx, dxi dx, dx, 


(A. 19) 


(A. 20) 
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where 


l> e, - p( + 

(A. 21 ) 

Favre weighting equation (A. 20) yields 


Dpti dpvjt, ()TTi p dujffij ()<i, 

()l 1 ()x r (h { (hi (hi 

(A. 22) 

where 


6/ = t + ^ Uk'Uk 

(A. 23) 

As was clone with the momentum equation we rewrite the filter of the product of two 

variables to obtain useable forms. 


ptpij = ptfUi t p{CtUj — t 1 u , ) 

(A. 24) 

pTij - fnii \ (W, - P u <) 

(A. 25) 

UjVjj = UjVij + ( U J a tj " U 3 (T ij) 

(A. 26) 

■£l 

II 

■£> 

+ 

1 

(A. 27) 

where 


- i 3t 

(A. 28) 

- i af 

« = - h ih. 

(A. 29) 
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and k = k(T). Substituting equations (A. 24) - (A. 27) into (A. 22) we obtain the 
following form of the filtered energy equation. 


Optt t ()pup, t c)u,p dll-fin Oq, () r 

+ — : + — — - —7—^ - 77— [p {( ,u- - c ,?/,)] 


i)t 


()-r> O.Vj Or, dx; dr, 


ft , . o . _ . d 

faT. ~ P M '-) + ^ (ujcr tj - UjtTij) (q; - q,) (A. 40) 


I he underbraced terms can be further manipulated to obtain terms for which sub-grid 
scale models have been previously developed. 

The argument of the derivative in term (?) is transformed as follows. 
p [t fU t — e,u t ) = p u t (c + ~upu k .) — u, (r -f ^upu k ) | 

= P [«iC + | ??,??*???• — Up + ^Uiujpi k ] (A. 31 ) 

= p{Ui'( — Up) + fp (uiUpUk — UjUpil k ) 

The first term in (A. 41) can be written as the sub-grid scale heat flux. 

p {up — up) = p (y U,c v T — UjC v T J 

= p£- {pr - u t f) (A. 32) 

= JLn. 

•>-1 

The second term in (A.31) is the sub-grid scale turbulent diffusion and is denoted as 

\p {uiUk.Uk - UiUpu. k ) = pD, 

Finally the argument of the derivative of term (?) in (A. 30) is simplified to 

R 


p{e,u, - e t Ui ) 

Using the filtered equation of state 


~Qi + pDj 


P = pRT 
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the argument of the derivative in term (») in ( A.30) can be rewritten as the sub-grid 


scale heat flux as was done in equation (A.d2). 

pu t — pu r = pRTiij — p H 7 Uj 


Term (m) in 


= pH (u,T - uj 7') 
= HQ, 

is transformed as follows. 


[A. 36) 


d , 

du ,o 

Oujdij 

[Uid 

dx , v J 

-■W)- dr . 

dx t 


du-jOij 

du-jOij 


dxj 

dx, 


do,, 

= u j .j 

dx, 

duj 

+ l7 -a.r, 


A. 31 


do 


— u , 




= e + u j - 


do. 


— a , 


1 dx, 
dd,j\ 


- O, 


dx, 


dr, J dx, 

where e is the sub-grid scale turbulent dissipation rate 


du, „ it; 

< -- °d- - a d — 


dx; dXi 

Substituting (A.34) (A. 36) and (A. 37) into (A.30) and rearranging terms gi 
the final form of the filtered energy equation. 


ves us 


dpt i dpuid duiP _ dtijBij dg, — d / 

— + “ + — - dXt Ox, dx, V-- 1 e '> dx 


dt dx i 


dx, 


dpDt 

dx, 


+ 


d ( do i 


n i ~Ts u :i «: 


do,, \ d „ 

- o— (?.■ - 9/) 


dx. 


dx, 


dx, V •' dx, 

A. 5 Determination of Pressure 

Pressure is normally obtained from the total energy as follows 

p = ph- l)(e, - \u k .u,) ( A.40) 
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Applying t he spatial filter and Favre weighting we obtain 

V = (? - 1 ) (W, ~ L 2 WkU k ) 


__ (A. 11 ) 

= b - J ) - ipukitk) 

r he above equation for pressure contains the filter ol a product of variables. iTpii/,. 
1 his undetermined quantity is eliminated by using the sub-grid scale kinetic energy 
p = b ~ 1 ) (pPi - l 2 Pmik) 

= (7 - 1 ) [pe, - ~pu k u k - \p (tij7u L - u k u k )] (A. 42) 

= (')'- 1 ) (ptt ~ \pPkUk - \t h .) 
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